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ABSTRACT :  Engineering  design  is  a  decision-making  process.  Optimization 
techniques  can  be  used  to  insure  that  better  decisions  are  made.  One  design  of  great 
interest  to  engineers  is  that  of  high-velocity  channels  used  for  routing  floodwater  out  of 
urban  areas.  In  the  design  of  these  channels  it  is  very  important  to  avoid  such  hydraulic 
phenomena  as  standing  waves,  hydraulic  jumps,  and  shocks.  These  will  require  higher 
wall  heights  and  more  expense.  These  channels  can  be  modeled  with  physical  models, 
but  they  are  expensive  and  time  consuming.  To  minimize  the  cost  of  building  and 
changing  the  physical  models  and  the  time  required  to  perform  the  study,  an  automated 
numerical  model  can  be  used  to  test  a  range  of  designs  before  construction  of  the  physical 
model.  The  resulting  design  can  be  used  as  an  initial  design,  which  is  close  to  the  desired 
design  requiring  fewer  changes  to  the  physical  model,  saving  time  and  money. 
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CHAPTER  I 


INTRODUCTION 

1.1  General 

The  population  of  Los  Angeles  County  had  grown  to  more  than  a  half  million 
residents  by  1910.  Few  had  experienced  the  prior  overflows  of  the  Los  Angeles 
River  and  even  fewer  understood  the  need  for  flood  control,  but  they  soon  would. 
On  February  18,  1914,  shortly  after  midnight,  rain  began  to  fall  on  the  already 
saturated  ground.  It  rained  for  three  days  delivering  more  than  19  inches  of  rain 
in  some  locations.  Los  Angeles  received  over  7  inches,  including  1.5  inches  in  an 
hour.  The  Los  Angeles  River,  the  San  Gabriel  River,  the  Rio  Hondo,  and  Ballona 
Creek  all  overflowed.  Flood  discharge  was  31,400  cubic  feet  per  second,  which 
is  equal  to  the  normal  flow  of  the  mighty  Colorado  River.  11,763  acres  in  Los 
Angeles  were  inundated  by  floodwaters.  The  flooding  destroyed  over  100  roads 
and  washed  out  35  bridges.  Rail  service  was  suspended  and  communication  with 
the  outside  world  was  interrupted  for  nearly  a  week.  Four  million  cubic  yards  of 
silt  were  dumped  into  Los  Angeles  and  Long  Beach  harbors  by  the  flood-swollen 
Los  Angeles  River,  rendering  some  channels  unnavigable.  Damage  throughout  the 
county  was  estimated  at  $10  million  ($151  million  in  1995  dollars).  Later,  county 
assessors  calculated  that  flooding  had  reduced  the  value  of  property  in  the  county 
by  $20  million.  The  devastating  flood  of  1914  clearly  demonstrated  that  the  Los 
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Angeles  River  could  no  longer  be  allowed  to  wander  at  will.  It  must  be  controlled. 
[Gumprecht  97] 

The  design  of  structures  to  control  waterways  such  as  the  Los  Angeles  River  is  a 
major  concern  for  engineers.  The  options  for  flood  control  in  urban  areas,  however, 
are  limited.  A  large  fraction  of  the  ground  surface  is  paved  causing  concentrated 
flood  flow  peaks.  It  is  not  feasible  to  build  a  reservoir  in  downtown  Los  Angeles  or 
modify  the  landscape  to  control  this  flow.  A  practical  method  of  routing  the  water 
through  the  urban  areas  is  via  the  use  of  high-velocity  channels.  High-velocity 
channels  offer  the  capability  to  carry  supercritical  flow. 

A  lot  of  details  must  be  considered  when  designing  flood  control  structures. 
Engineering  design  is  a  decision-making  process.  Optimization  techniques  can 
aid  in  ensuring  that  better  decisions  are  made  concerning  the  design.  An  area 
of  engineering  design  that  can  benefit  from  optimization  techniques  is  the  design 
and  modification  of  high-velocity  channels  essential  for  the  routing  of  floodwater 
through  urban  areas.  The  proper  design  of  new  channels  and  re-design  of  existing 
channels  is  required  to  avoid  such  things  as  bank  erosion,  damaged  equipment, 
increased  operating  expenses,  flooding,  and  higher  construction  costs.  Physical 
models  are  very  useful  as  a  tool  to  determine  appropriate  designs  of  flood  channels 
to  meet  certain  site-specific  criteria,  but  initial  design  and  modifications  of  physical 
models  are  very  costly  in  both  time  and  money.  Due  to  the  time  and  cost 
constraints  of  physical  models,  it  is  not  practical  to  examine  a  wide  range  of  designs. 
This  could  result  in  hydraulic  performance  that  is  only  acceptable  over  a  limited 
range.  The  ideal  scenario  is  to  build  the  initial  physical  model  as  close  to  the  final 
design  as  possible.  This  can  be  accomplished  by  using  an  automated  hydraulic 
design  program  to  produce  a  better  design  prior  to  the  building  of  the  physical 
model  which  will  reduce  both  design  time  and  cost. 
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Burg  [Burg  99]  has  built  such  a  capability.  He  explored  numerous  optimization 
methods  to  determine  those  most  suitable  for  problems  involving  high-velocity 
channels.  Burg’s  work  demonstrates  an  optimization  strategy  integrated  with 
the  HIVEL2D  analysis  code  for  design  optimization  of  simple  two-dimensional 
test  cases.  HIVEL2D,  a  two-dimensional  numerical  model  that  solves  the  shallow 
water  equations,  is  a  proven  tool  for  modeling  these  types  of  problems  [Berger  95]. 
Examples  presented  by  Burg,  however,  were  represented  by  simple  geometries.  The 
examples  deal  with  flume  type  problems,  as  well  as  super  elevation  in  circular  bends 
and  embedded  bodies  such  as  bridge  piers.  The  application  of  this  technique  to  an 
urban  flood-control  channel,  previously  designed  through  an  experimental  study, 
is  the  concern  of  this  investigation.  Application  to  such  a  “real-world”  problem 
will  provide  an  assessment  and  demonstration  of  the  utility  of  the  design  method. 

1.2  Objectives 

The  objective  of  this  work  is  to  assess  the  practicality  of  using  the  optimization 
technique  developed  by  Burg  to  aid  in  the  design  of  a  realistic  high-velocity  channel. 
Since  the  test  cases  used  by  Burg  were  represented  by  simple  geometries  and  did 
not  address  the  more  complex  features  of  channels,  the  steps  required  to  address 
more  complex  problems  must  be  established.  Therefore,  the  procedure  for  applying 
design  optimization  techniques  to  “real-world”  problems  will  be  determined  by  this 
research.  The  desired  outcome  of  the  research  is  to  apply  Burg’s  technique  to  a 
“real-world”  problem  and  to  develop  a  procedure  for  similar  future  applications. 
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1.3  Scope 

A  numerical  model  of  the  Walnut  Creek  physical  model  is  developed  and 
simulation  results  are  compared  to  the  physical  model  results.  A  series  of  model 
parameters  are  tested  to  determine  the  model  sensitivities.  This  reduces  the 
number  of  parameters  to  only  those  that  have  a  major  impact  on  the  design. 
The  results,  along  with  engineering  judgment,  are  used  to  explore  and  determine 
the  appropriate  parameters  to  be  optimized  and  to  explore  the  applicability  of 
possible  objective  functions  to  be  used  in  the  model.  Once  the  appropriate  design 
parameters  and  objective  function  are  determined,  the  design  is  automated  to 
yield  a  better  hydraulic  design  for  the  given  constraints.  Walnut  Creek  was  used 
as  an  evaluation  case  since  the  design  of  this  hydraulic  structure  was  extensively 
analyzed  via  an  experimental  project.  [Davis  87] 


CHAPTER  II 


BACKGROUND 

2.1  Modeling  High-Velocity  Channels 

High-velocity  channels  are  used  for  drainage  in  urban  regions,  since  urban 
sprawl  increases  rainfall  runoff  due  to  altered  land  use.  Flood  control  channels 
are  designed  and  built  to  safely  manage  the  anticipated  hydrologic  load.  The 
desire  is  to  minimize  the  water’s  time  of  residence  in  the  urban  area.  The  channels 
are  designed  to  carry  supercritical  flow  to  reduce  the  water  depths  and  the  required 
route.  Structures,  such  as  bends,  transitions,  and  confluences  cause  flow  to  choke, 
form  jumps,  and/or  form  standing  waves.  These  hydraulic  conditions  generally 
necessitate  higher  walls,  bridges,  and  other  costly  containment  structures.  A 
poorly  designed  channel  can  cause  bank  erosion,  damaged  equipment,  increased 
operating  expense,  and  reduced  efficiency  [Berger  95].  Furthermore,  crossings  may 
be  washed  out,  and  the  town  may  flood.  To  improve  the  design  of  channels, 
engineers  use  models  to  reproduce  the  channels  and  run  a  suite  of  tests  to  determine 
the  functionality  of  the  channels  under  selected  hydraulic  conditions.  The  two  ways 
of  modeling  these  channels  are  with  physical  models  and  with  numerical  models. 

2.1.1  Physical  Models 

Numerical  modeling  of  high-velocity  channels  has  emerged  as  a  tool  to  augment 
physical  models.  Physical  models  are  not  as  popular  as  they  once  were,  but  they  are 

still  very  useful  in  channel  design  due  to  the  limitations  found  in  numerical  models. 
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One  such  limitation  of  shallow  water  models  is  that  they  are  typically  hydrostatic. 
Physical  models  can  model  things  such  as  non-hydrostatic  flow  at  the  tail  of  a 
bridge  pier  that  a  shallow-water  numerical  model  is  not  capable  of  modeling. 
Physical  models  are  extremely  expensive  to  build,  costing  on  the  order  of  $150,000 
per  study  and  are  time  consuming  requiring  an  average  of  three  months  to  complete 
[Stockstill  00].  Changes  to  the  physical  model  require  a  “cut  and  try”  technique 
that  involves  tearing  down  the  unwanted  sections  of  the  channel  and  rebuilding 
them  with  the  new  desired  design  [Soulis  92],  Each  modification  on  average  costs 
$30,000  and  take  a  minimum  of  1  month  to  complete.  Though  physical  models  can 
reproduce  details  of  actual  hydraulic  structures,  they  are  subject  to  the  limitation 
of  scale  modeling.  It  is  impossible  to  reproduce  the  physical  problem  to  scale. 
The  problem  is  usually  scaled  down  to  a  more  cost  efficient  and  manageable  size. 
For  instance,  a  physical  model  constructed  on  a  scale  of  1:25  means  that  1  foot 
of  model  is  equal  to  25  feet  of  prototype.  Scaling  the  model  to  a  more  reasonable 
size  reduces  the  space  and  materials  required  to  build  the  physical  model.  The 
equations  of  hydraulic  similitude,  based  on  Froudian  relations,  are  used  to  express 
mathematical  relations  between  the  dimensions  and  hydraulic  quantities  of  the 
model  and  prototype.  The  bed  slope  of  the  physical  model  is  adjusted  to  accoimt 
for  the  difference  in  the  roughness  of  the  model  material  and  the  prototype.  Table 
2.1  gives  the  scale  relation  between  the  physical  model  and  prototype. 

A  physical  model  study  conducted  at  the  Waterways  Experiment  Station  of 
the  Walnut  Creek  Flood-Control  Project  in  Contra  Costa  County,  California,  is 
considered  for  testing  the  applicability  of  Burg’s  techniques  to  an  actual  problem. 
The  physical  model  study  was  a  1:25  reproduction  of  approximately  1,084  ft  of 
the  San  Ramon  Bypass  Channel,  730  ft  of  the  Walnut  Creek  Channel  upstream 
from  the  jimction,  and  640  ft  of  the  Walnut  Creek  Channel  downstream  from 


Table  2.1:  Scale  Relations  Between  the  Physical  Model  and  Prototype 


Characteristic 

Dimension  (terms  of  length) 

Model:  Prototype 

Length 

J_jtp  —  Xj 

1:25 

Area 

A  =  L2r 

1:625 

Velocity 

vr  =  l) 

1:5 

T  ime 

Tr  =  Lj 

1:5 

Discharge 

Qr  =  Lr 

1:3,125 

Roughness  Coefficient 

Nr  =  L} 

1:1.71 

the  junction  (Figure  2.1).  The  project  was  aimed  at  making  improvements  to 
the  channel  to  provide  flood  protection  to  about  6,670  acres  in  the  floodplain  at 
and  below  the  city  of  Walnut  Creek.  The  site  contains  several  hydraulic  structures 
including  circular  curves  with  super  elevation  and  transition  spirals,  a  divider  wall, 
a  contraction,  and  a  confluence.  Figure  2.3  shows  the  layout  of  the  channels  and 
Figure  2.4  shows  a  general  view  of  the  physical  model. 

A  series  of  floods  had  necessitated  re-evaluation  of  the  Walnut  Creek  Flood- 
Control  channel  design.  In  1955,  Walnut  Creek  had  a  peak  discharge  of  11,000 
cfs  and  San  Ramon  Creek  had  a  peak  discharge  of  6,900  cfs.  As  a  result,  3,500 
acres  flooded  and  1,000  homes  and  50  businesses  were  affected  with  the  damage 
estimated  to  be  $1  million.  A  flood  in  1958  had  similar  results.  The  peak  discharge 
of  Walnut  Creek  was  12,200  cfs.  3,400  acres  were  flooded  and  140  businesses  and 
800  homes  were  damaged  as  shown  in  Figure  2.2.  The  estimated  cost  of  damages 
was  $l|  million.  In  1960,  the  Flood  Control  Act  was  passed  which  authorized 
the  Walnut  Creek  Basin  Project.  Smaller  floods  followed  in  1962  and  1963.  The 
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100-year  event  was  estimated  in  1964  to  be  18,000  cfs.  Design  improvement  for 
the  Walnut  Creek  Basin  soon  began.  In  the  late  1970s,  however,  flood  problems 
continued  in  the  San  Ramon,  Las  Trampas,  and  Walnut  Creeks  due  to  the  small 
capacity  of  the  existing  channels.  Urbanization  had  decreased  infiltration  and 
increased  runoff  in  the  area.  In  1982,  peak  flow  of  7,000  cfs  occurred  in  San  Ramon 
Creek  causing  flow  to  approach  the  top  of  the  creek  banks  and  lap  at  bridges.  The 
solution  was  determined  to  be  a  bypass  channel  connecting  San  Ramon  Creek  and 
Walnut  Creek.  The  purpose  of  the  physical  model  study  was  to  determine  the 
adequacy  of  the  channel  during  100-year  frequency  flow  conditions  and  to  develop 
modifications  to  increase  the  hydraulic  capacity  and  improve  flow  conditions.  That 
is,  minimize  cross-waves  generated  by  the  curves  and  the  jimction. 

The  main  focus  of  the  study  was  the  length  of  the  divider  wall  extension  at  the 
junction.  Channel  junctions  are  one  of  the  more  important  hydraulic  problems, 
because  as  flows  from  smaller  channels  are  combined  standing  waves  may  be 
produced  necessitating  increased  wall  heights  in  the  vicinity  of  the  junction.  Also, 
hydraulic  jumps  may  form  in  one  or  both  of  the  channels  if  the  divider  wall  is  not 
correctly  located  or  scaled.  Figure  2.5  shows  a  hydraulic  jump  and  standing  waves 
in  the  vicinity  of  the  junction.  Several  divider  wall  extension  lengths  were  modeled 
under  various  flow  conditions.  Testing  of  divider  wall  extensions  are  simpler  than 
testing  other  parameters  because  it  only  requires  that  sections  of  different  lengths 
be  placed  in  the  model  and  oriented  by  the  engineer.  There  is  no  need  for  tearing 
down  and  rebuilding  the  model.  Other  areas  of  interest  were  the  width  of  the  San 
Ramon  Bypass  Channel  and  improved  channel  alignment  at  the  junction. 

In  1984,  the  estimated  100-year  flood  event  for  the  Walnut  Creek  Flood-Control 
Project  was  22,000  cubic  feet  per  second.  Currently,  the  100-year  flood  event  has 
not  changed  from  its  estimated  value  in  1984.  Because  the  San  Ramon  Bypass 
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Figure  2.1:  Vicinity  map  of  the  Walnut  Creek  Physical  Model  Study. 
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Figure  2.2:  Flood  of  1958  in  Walnut  Creek,  California 


Channel  connects  San  Ramon  Creek  and  Walnut  Creek  arid  bypasses  the  junction 
between  Las  Trampas  Creek  and  San  Ramon  Creek  which  forms  Walnut  Creek, 
the  amount  of  flow  in  the  San  Ramon  Bypass  Chaimel  effects  the  amount  of  flow 
in  Walnut  Creek.  Two  discharge  conditions  were  considered.  The  first  condition 
allowed  for  maximum  flow  in  the  San  Ramon  Bypass  Channel  with  the  concurrent 
flow  set  for  Walnut  Creek.  The  second  condition  allowed  for  the  maximum  flow  in 
Walnut  Creek  with  the  concurrent  flow  set  for  the  San  Ramon  Bypass  Chaimel.  In 
prototype,  increasing  discharges  caused  standing  waves  to  develop  downstream  of 
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Figure  2.4:  General  View  of  the  Physical  Model  [Davis  87]. 
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Figure  2.5:  Standing  Waves  and  Hydraulic  Jump  near  Junction  [Davis  87]. 
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the  confluence  and  extend  several  hundred  feet  due  to  the  poor  channel  alignment 
as  illustrated  in  Figure  2.6.  Water  depth  exceeded  the  wall  heights  at  several 
points  downstream  from  the  junction  as  illustrated  in  Figure  2.7.  Adding  a  divider 
wall  extension  caused  flow  to  choke  and  form  a  hydraulic  jump  in  Walnut  Creek 
just  before  the  junction  that  spread  upstream  as  the  length  of  the  divider  wall 
was  increased.  Figure  2.8  shows  the  jump  pushed  upstream  in  Walnut  Creek 
about  70  feet  from  the  junction.  The  longer  the  divider  wall  extension,  however, 
the  smoother  the  water  surface  downstream  from  the  junction.  The  physical 
model  study  concluded  that  a  40-foot  long  divider  wall  extension  should  be  added, 
the  width  of  San  Ramon  Bypass  Channel  should  be  decreased,  and  the  channel 
alignment  at  the  junction  should  be  improved  [Davis  87]. 

The  number  of  designs  that  can  be  explored  by  a  physical  model  is  limited  by 
the  construction  time  and  cost,  though  defining  the  details  of  the  model  requires 
exploration  of  a  vast  range  of  alternatives.  A  numerical  model  such  as  HIVEL2D 
can  be  used  to  assess  the  design  computationally  before  construction  of  the  physical 
model  begins  and  to  screen  alternatives.  Using  an  automated  hydraulic  design 
package  would  accelerate  this  screening  process  and  lead  to  an  improved  initial 
physical  model  thus  reducing  the  time  spent  on  the  physical  model.  This  would 
allow  for  exploration  of  more  design  alternatives  in  a  shorter  length  of  time  resulting 
in  a  more  cost-effective  solution.  Though  the  HIVEL2D  numerical  model  cannot 
show  such  hydraulic  phenomena  as  non- hydrostatic  flow  at  the  tail  of  a  bridge  pier, 
which  would  require  a  physical  model  or  more  complex  computational  model,  it  can 
give  an  improved  initial  design  and  reduce  the  number  of  modifications  required. 
This  is  a  restriction  on  the  hydraulics  only  not  the  optimization. 
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Figure  2.6:  Standing  Waves  and  Overtopping  Downstream  of  Junction  [Davis  87]. 
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Figure  2.7:  Depths  Exceeding  Wall  Heights  Downstream  of  Junction  [Davis  87]. 
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Figure  2.8:  Hydraulic  Jump  Pushed  Out  of  View  [Davis  87]. 
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2.1.2  Numerical  Models 

Numerical  simulations  are  advantageous  because  they  are  relatively  fast  and 
inexpensive,  which  allows  for  more  changes  and  testing  of  more  parameters 
compared  to  physical  modeling.  Since  the  vertical  accelerations  are  usually  small, 
the  primary  basis  for  simulating  flow  in  high-velocity  channels  is  the  shallow 
water  equations.  The  shallow  water  equation  models,  however,  are  limited  by 
the  assumptions  made  in  the  equations  and  by  the  limitations  of  the  computer. 
The  size  of  the  computer  can  limit  the  amount  of  spatial  resolution  that  can  be 
used  to  represent  the  problem. 

There  are  numerous  numerical  models  for  solving  open- channel  flow.  For 
example,  Zhou  and  Stansby  [Zhou  99]  developed  a  model  to  predict  the  hydraulic 
jump  in  a  straight  open  channel.  Guillou  and  Nguyen  [Guillou  99]  present  a 
technique  for  solving  the  two-dimensional  shallow  water  problems  based  on  a  finite 
volume  technique.  The  model  chosen  for  this  study  is  a  finite  element  model, 
HIVEL2D,  developed  by  Berger  and  Stockstill  [Berger  95]. 

HIVEL2D  is  appropriate  for  this  work  because  of  its  capability  to  model 
supercritical,  subcritical,  and  transition  flow.  Since  the  objective  of  this  study  is  to 
identify  non-smooth  flow  or  sloshing,  a  model  of  at  least  two-dimensions  is  required. 
The  vertical  accelerations  for  this  problem  are  negligible  compared  to  the  effects  of 
gravity.  HIVEL2D  solves  the  shallow  water  equations  in  conservative  form,  which 
conserves  momentum  and  allows  for  the  correct  calculation  of  the  shock  speed  and 
location.  Three-dimensional  models  can  be  used  for  this  simulation,  but  they  are 
more  costly  and  have  capabilities  that  are  not  required  to  solve  this  problem.  The 
bed  slope  in  the  Walnut  Creek  problem  is  geometrically  mild  though  hydraulically 
steep,  which  meets  the  assumption  made  by  HIVEL2D.  HIVEL2D  is  designed  for 
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solving  problems  with  rapidly  varying  flow  and  is  appropriate  for  the  Walnut  Creek 
problem. 


2.2  Literature  Review 

Since  optimization  can  be  a  powerful  tool  for  engineers  during  the  design 
process,  optimization  techniques  are  used  in  a  vast  range  of  engineering  areas. 
Numerous  optimization  techniques  have  been  applied  to  engineering  problems  in 
the  areas  of  aerospace,  groundwater,  and  surface  water. 

2.2.1  Aerospace  Applications 

Makinen  et  ah  [Makinen  99]  applied  optimization  to  a  multi-objective, 
multidisciplinary  design  optimization  of  a  two-dimensional  airfoil.  Genetic 
algorithms  were  used  to  obtain  an  approximation  for  the  Pareto  set  of  optimal 
solutions.  This  work  is  not  traditional  because  it  optimizes  more  than  one 
objective  function  of  more  than  one  discipline.  The  two  objective  functions  are 
the  drag  coefficient  with  constraints  of  the  lift  coefficient  above  a  given  value  with 
CFD  analysis  solvers  based  on  finite  volume  discretizations  of  the  inviscid  Euler 
equations  and  the  integral  of  the  transverse  magnetic  radar  cross  section  over  a 
given  sector.  Several  non-dominated  designs  were  obtained. 

Design  sensitivity  analysis  is  discussed  in  Hou  et  al.  work  [Hou  94],  This 
method  is  applied  to  aerodynamic  problems.  A  derivation  of  the  aerodynamic 
sensitivity  equations  is  presented  that  uses  the  existing  formats  of  implicit 
algorithms  for  solving  the  Euler  equations.  The  derived  sensitivity  equations  are 
efficiently  solved  with  a  direct  equation  solver  for  small-scale  problems. 

Optimization  can  be  used  to  determine  airfoil  design  such  that  pre-specified 
design  criteria  are  met.  Soemarwoto  et  al.  [Soemarwoto  99]  explored  the  area  of 
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aerodynamic  airfoil  shape  optimization.  This  work  applied  the  adjoint  operator 
approach,  utilizing  a  compressible  inviscid  flow  model  based  on  the  Euler  equations 
and  a  compressible  viscous  flow  model  based  on  Reynolds-averaged  Navier  Stokes 
equations.  The  model  solves  constrained  transonic  aerodynamic  design  (pressure 
drag  reduction)  problems. 

The  emphasis  on  reduction  of  design  cycles,  time,  and  cost  in  design  of 
commercial  aircraft  has  sparked  a  renewed  interest  in  design  optimization  in 
aerodynamic  structures  and  aeroelastics.  In  their  work,  Melvin  et  al.  [Melvin  99] 
use  the  TRANAIR  code  to  optimize  design.  TRANAIR  is  a  two-  and  three- 
dimensional  full  potential  flow  code  with  directly  coupled  strip  boundary  layer  and 
is  capable  of  handling  complex  geometries  through  solution  adaptive  local  grid 
refinement.  The  sensitivity  method  is  preferred  to  the  adjoint  method  because 
some  second-order  information  is  naturally  and  inexpensively  available  in  the 
sensitivity  method.  The  problem  is  solved  using  Newton’s  method. 

2.2.2  Groundwater  Applications 

Optimization  methods  are  also  very  useful  for  solving  groundwater  problems. 
Townley  et  al.  [Townley  85]  apply  present  computational  algorithms  to  steady 
and  transient  models  for  groundwater  flow.  The  aquifer  storage  coefficients, 
transmissivities,  distributed  inputs,  and  boundary  values  may  all  be  uncertain. 
The  discrete  derivation  of  the  adjoint  method  is  used  to  calculate  the  derivative  of 
a  function  with  respect  to  the  parameters  of  the  transient  numerical  flow  models. 
The  explicit  calculation  of  sensitivities  and  the  calculations  of  the  gradient  of  the 
objective  function  in  the  inverse  problem  are  calculated  with  the  adjoint  method. 
A  Gauss- Newton  line  search  procedure  is  used  to  relate  the  second  derivative  in  the 
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direction  of  a  search  to  the  projection  of  the  sensitivity  matrix  onto  the  direction 
of  search. 

Adjoint  state  equations  are  used  in  Sun  et  al.  [Sun  92]  work  for  the  stochastic 
partial  differential  equations  relating  transient  head  and  log  hydraulic  conductivity 
perturbations.  The  model  reliability  is  evaluated  through  the  variance  estimate 
method  using  adjoint  sensitivity  analysis  and  the  cokriging  estimate.  A  stochastic 
inverse  procedure  for  transient  groundwater  flow  was  developed.  By  using  the 
adjoint  state  equation  and  solving  for  the  expected  head,  all  elements  of  the 
covariance  matrix  can  be  calculated  based  on  the  first-order  approximation. 

Guan  et  al.  [Guan  99]  use  optimization  to  design  a  pump-and- treat 
groundwater  remediation  system.  A  new  computational  procedure  called  a 
progressive  genetic  algorithm  is  used  to  minimize  the  total  cost  of  the  pump- 
and-treat  system  while  defining  the  locations  and  extraction  or  injection  rate  of 
the  wells.  Constraints  of  concentration,  velocity,  and  injection  and  extraction 
balance  were  considered.  These  techniques  proved  reliable  and  robust  for  the 
given  problem. 

Optimization  techniques  have  also  been  used  to  locate  wells  for  monitoring 
purposes.  In  work  done  by  Storck  et  al.  [Storck  97],  the  flow  and  transport 
simulations  are  passed  to  an  optimization  model.  The  optimization  problem  is 
solved  by  simulated  annealing.  The  problems  consist  of  three  conflicting  objectives; 
maximum  detection  probability,  minimum  cost  (number  of  wells),  and  minimum 
volume  of  contaminated  groundwater  at  the  time  of  detection.  Application  of  this 
model  to  groundwater  problems  show  that  the  trade-off  curves  for  the  objectives 
based  on  too  few  realizations  over  predicts  monitoring  performance. 
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2.2.3  Surface  Water  Applications 

Atanov  et  al.  [Atanov  98]  developed  an  optimization  method  for  minimizing 
water-level  fluctuations  in  an  open-channel  controlled  by  pumping  stations  on  both 
ends.  The  problem  was  a  variational  problem  to  determine  the  optimum  flow 
control,  given  constraints  at  the  opposite  end.  A  global  measure  of  the  water  level 
deviation  away  from  a  desired  water  level  integrated  in  space  and  time  was  used 
as  the  objective  function.  The  goal  was  to  minimize  deviations  while  satisfying 
governing  equations  and  boundary  conditions.  The  problem  was  solved  using  a 
direct  flow  solver  and  a  conjugate  determination  of  the  Lagragian  multiplier  solved 
backward  in  time. 

Hsu  et  al.  [HsuM  99]  used  a  vertical  (laterally  averaged)  two-dimensional 
model  of  a  branched  estuarine  river  system  to  determine  friction  and  turbulent 
diffusion/dispersion  coefficients.  The  coefficients  affect  the  calculations  of  the 
surface  elevation  velocity  and  salinity  distribution. 

Soulis  et  al.  [Soulis  92]  designed  open-channel  expansions  and  contractions  to 
give  specified  distributions  of  depth.  Prescribing  depth- velocity  values  can  avoid 
boundary  layer  separation  and  cavitation.  The  usual  design  process  involves  the 
use  of  physical  and/or  numerical  models  to  resolve  the  fluid  properties  adequately. 
The  combined  effects  of  physical  and  numerical  approaches  can  lead  to  designs 
acceptable  from  an  engineering  point  of  view.  Computer-aided  design  will  improve 
the  performance  of  any  hydraulic  structure  while  shortening  the  amount  of  human 
effort  needed  for  such  a  design.  The  authors  developed  a  general  numerical  method 
of  design  for  two-dimensional  channel  expansions  and  contractions  with  prescribed 
depths  and  velocities  along  the  channel  walls  using  a  finite  volume  analysis  code. 
The  process  iterates  between  the  direct  solution  and  an  inverse  solution.  For 
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a  typical  problem  the  procedure  should  be  incorporated  20-30  times  in  order  to 
converge  to  an  acceptable  geometry.  This  procedure  essentially  reverses  the  roles  of 
the  independent  variables  x  and  y  (spatial  variables)  and  the  dependent  variables 
h,  p,  and  q  (flow  variables)  at  specific  locations.  The  application  of  this  method 
is  limited  and  not  as  general  as  the  method  used  in  this  work. 

Mohammadi  et  al.  [Mohammadi  99]  used  an  optimal  shape  design  and 
unstructured  mesh  deformations,  automatic  differentiation  for  the  gradient 
computation,  and  mesh  adaption  by  metric  control  in  two-dimensions.  For  a  CAD- 
free  design,  the  only  geometrical  entity  available  during  optimization  should  be 
meshed.  The  model  uses  the  reverse  mode  of  automated  differentiation  to  produce 
gradients  of  discrete  operators.  The  mesh  is  part  of  the  optimization  procedure 
and  a  function  of  the  solution.  This  avoids  mesh  dependencies  in  optimization  and 
the  direct  problem. 

In  work  done  by  Zhu  et  al.  [Zhu  99]  optimization  techniques  were  used  to 
determine  the  optimal  locations  and  scheduling  of  dredging  to  minimize  cost  and 
obtain  a  channel  where  water  depths  are  not  less  than  some  specified  depth.  The 
control  of  sedimentation  can  be  extremely  costly;  hence  optimization  is  a  useful 
tool.  The  optimization  problem  is  solved  by  the  conjugate  gradient  method.  The 
gradient  of  the  cost  function  (objective  function)  is  calculated  by  solving  the  adjoint 
problem.  Through  this  study  it  was  determined  that  for  a  nonlinear  model  the 
adjoint  approach  is  only  valid  for  small  variations  in  dredging  depth. 

Mousavi  and  Ramamurthy  [Mousavi  00]  introduced  a  new  composite  algorithm 
for  optimizing  the  operating  policy  of  multi- reservoir  systems.  The  model  uses 
optimal  control  theory  and  penalty  successive  linear  programming  as  techniques 
for  modeling  large  and  complex  water  supply  systems.  The  objective  function 
minimizes  the  required  reservoir  costs  to  supply  specified  yields. 
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In  their  work,  Cunha  et  al.  [Cunha  99]  use  the  Newton  search  method  to  solve 
the  hydraulic  network  equations  to  obtain  the  least-cost  design  of  a  looped  water 
distribution  network  using  simulated  annealing.  Simulated  annealing  works  well 
for  large-scale  optimization  problems  that  are  cast  in  discrete  or  combinatorial 
form.  The  optimal  water  distribution  design  problems  discussed  in  the  study 
contain  discrete  elements.  To  keep  these  designs  close  to  reality  they  must  be 
formulated  as  large  size,  nonlinear  mixed  integer  models.  Results  show  that  this 
method  can  provide  high  quality  solutions  for  network  design  problems  compared 
to  past  studies. 

Piasecki  et  al.  [Piasecki  97]  used  adjoint  sensitivity  analysis  to  determine  the 
control  of  contaminant  releases  in  rivers.  Two-dimensional  finite-element  models 
are  used  to  determine  the  hydrodynamic  and  mass  transport  conditions.  The 
sensitivities  of  the  loading  parameters  are  computed.  This  method  reproduces 
almost  identical  “sensitivity  functions”  with  repeated  numerical  solutions. 

River  and  channel  flood  routing  are  modeled  using  the  nonlinear  Muskingum 
models  in  Mohan’s  work  [Mohan  97].  Mohan  used  genetic  algorithms  to  estimate 
the  parameters  of  the  Muskingum  models.  Comparisons  of  the  performance  of  the 
genetic  algorithm  and  other  parameter  estimation  procedures  were  evaluated.  The 
genetic  algorithm  estimates  the  parameters  quickly  and  objectively  and  performed 
better  or  as  well  as  other  methods  to  which  it  was  compared. 

Sanchez  et  al.  [Sanchez  98]  discuss  the  use  of  neural  networks  as  a  means  of 
reducing  computing  costs  of  coastal  sewage  systems.  Neural  network  calculations 
are  used  to  reduce  calculation  time  for  studying  the  influence  of  storm  discharge 
on  the  bacteriological  quality  of  bathing  waters.  Neural  networks  maintain  an 
approximation  level  similar  to  numerical  solutions.  The  study  concludes  that 
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neural  network  with  optimized  functional  link  is  suitable  for  design  of  the  collector 
network  of  coastal  sewage  systems  when  storm  discharges  occur. 

Burg  [Burg  99]  developed  an  optimization  strategy  for  open-channel  flow 
problems.  In  his  work,  optimization  methods  were  analyzed  to  determine  their 
applicability  to  various  problems.  Table  2.2  shows  a  comparison  of  the  various 
optimization  methods.  The  code  used  in  this  work,  HIVEL2D,  is  a  high-fidelity 
implicit  code.  From  the  table,  the  discrete  sensitivity  analysis  is  reported  to  be 
best  for  this  type  of  problem,  whereas  the  continuous  sensitivity  analysis  is  suitable 
for  explicit  high-fidelity  problems.  The  inverse  method  and  genetic  algorithms  are 
not  appropriate  for  this  type  of  problem.  The  adjoint  variable  formulation  or  the 
quasi-analytic  method,  described  below,  may  be  used  to  compute  the  design  space 
gradients.  Updating  the  design  variable  is  accomplished  by  the  method  of  steepest 
descent  with  a  step  size,  the  method  of  steepest  descent  with  a  linear  search,  the 
Gauss-Newton  algorithm,  or  the  BFGS  method.  The  technique  was  applied  to 
channel  expansion,  channel  contraction,  and  embedded  body  problems. 

The  adjoint  method  was  used  to  calculate  the  design  space  derivative  of  the 
objective  function,  F(Q((3),  X(P),  P)  where  Q  is  the  steady-state  flow  variables,  X 
is  the  discretized  grid,  and  (5  is  the  vector  of  design  variables  or 

dF  _  dFT  dQ  dFT  dX  OF 
d(3i  dQ  dPi  +  dX  d/3i  +  d Pi 

In  this  equation  js  the  total  variation  of  F  with  respect  to  the  design  variable 
Pi.  |||,  and  can  be  calculated  since  the  dependency  of  F  on  Q,  x,  and  P  is 
explicit.  can  be  estimated  by  finite  differencing  the  results  of  the  grid  generation 
code  or  by  differentiating  the  explicit  dependencies  between  the  grid  and  the  design 
variables  with  hand- differentiation,  complex  Taylor  series  expansion,  or  ADIFOR 


27 


Table  2.2:  Comparison  of  Optimization  Methods  [Burg  99] 


Optimization 

Method 

Type  of 
Problem 

Advantages 

Disadvantages 

Inverse 

Methods 

Analytic 

Formula 

Highly 

Efficient 

Not  Generally 
Applicable 

Genetic 

Algorithm 

(Probabilistic 

Methods) 

Discontinuous, 

Discrete, 

Cheap  Simulations, 
Multi-Modal 

Avoids  Local 
Minima, 

No  Gradient 
Needed 

Many 

Function 

Evaluations 

Finite 

Difference 

Any 

Easiest 

To  Use 

Large  Computer 
Cost,  Accuracy 

ADIFOR, 

CTSE 

Any 

Highly  Accurate 
Derivative, 
Easy  to  Use 

Large 

Computational 

Cost 

Continuous 

Sensitivity 

Analysis 

Explicit 

High-Fidelity 

Computationally 

Efficient 

Derive  and 
Solve 

Adjoint  Equations 

Discrete 

Sensitivity 

Analysis 

Implicit 

High-Fidelity 

Accurate 

Derivatives, 

Efficient 

Jacobian 

Matrix 

Needed 

[Bischof  92],  Estimation  of  the  vector  ^  with  finite  differences  would  require 
an  additional  steady-state  simulation.  This  term,  however,  also  appears  in  the 
derivative  of  the  discretized  system  of  governing  equations,  W(Q(/3),  x(/3),  (3)  =  0 
or 


dW  _  dVVdQ_  m^d\_  dW_ 
df3i  d Q  d(3i  +  d\  dfa  +  <9/3* 


0 


(2.2) 


As  in  equation  2.1,  all  terms  except  for  can  be  calculated  without  the  need  for 
a  steady-state  simulation.  Equation  2.2  and  Equation  2.1  are  used  to  estimate  the 
design  space  derivative 

For  the  adjoint  variable  formulation  of  discrete  sensitivity  analysis,  equation 
(2.2)  is  multiplied  by  an  adjoint  vector  A  and  added  to  equation  (2.1).  This  yields 


(IF_  _  fdF^  TdlE\  o >Q  f dF^  TdVE\  dF_  TdW 
d(3i  \  dQ  +  d Q  )  d(3i+  \  d\  +  dx)  Ofa  +  <9$  +  dfa 


To  avoid  having  to  solve  for  the  coefficient  of  is  set  equal  to  zero.  That  is 


dFT 

~dQ 

dWT 

~w 


+ 

A 


AJ 


■  dW 
~dQ 
OF 
~dQ 


=  0 


(2.4) 


A  must  only  be  calculated  once  for  each  objective  function  and  is  independent 
of  the  design  variables  for  the  adjoint  variable  formulation.  The  adjoint  variable 
formulation  scales  with  the  number  of  functions  and  constraints  and  is  appropriate 
for  problems  with  only  one  objective  function  and  a  few  constraints.  The  direct 
formulation  scales  with  the  number  of  design  variables. 

For  the  direct  formulation  of  discrete  sensitivity  analysis,  or  quasi-analytic 
method,  the  vector  is  obtained  by  solving  equation  (2.2)  directly  or 

dW  dQ  _  _dW  dX  _  dW  _  _dW 

dQ  d(3i  dx  df3i  d(3i  df3i 

The  vector  ^  must  be  calculated  for  each  design  variable.  This  makes  the 
computational  cost  of  the  direct  formulation  greater  than  that  of  the  adjoint 
variable  formulation  when  the  number  of  design  variables  is  greater  than 


O  fixed 


(2.5) 
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the  number  of  objective  function  and  constraints.  However,  when  using  the 
Gauss-Newton  optimization  algorithm  the  number  of  objective  functions  to  be 
differentiated  is  much  larger  than  the  number  of  design  variables  since  this 
algorithm  views  each  term  in  the  objective  function  as  a  separate  residual  function. 
For  this  reason,  the  quasi-analytic  method  is  more  cost  effective.  [Burg  99] 

The  work  presented  in  this  research  deals  with  one  and  two  design  variable 
problems.  The  objective  function  is  in  the  non-linear  least  squares  form,  so 
the  Gauss-Newton  method  may  be  used.  The  quasi-analytic  method  is  used  to 
calculate  the  design  space  derivative. 


CHAPTER  III 


METHODOLOGY 

Burg’s  work  was  applied  to  several  types  of  high-velocity  channel  transitions. 
He  did  not  address  issues  such  as  the  confluence  of  two  channels,  channels  with 
multiple  design  flows,  or  channels  with  multiple  transition.  The  purpose  of  this 
research  is  to  expand  Bnrg’s  tests  to  include  those  issues  not  addressed  by  Bnrg 
by  applying  his  work  to  a  “real-world”  problem.  The  Walnut  Creek  Flood-Control 
Channels  are  high-velocity  channels  that  include  hydraulic  structures  such  as 
contractions,  bends,  and  a  confluence.  The  design  flows  from  the  two  channels 
differ  by  as  much  as  8,400  cfs,  and  the  channels  contain  both  sub-  and  supercritical 
flows  with  a  hydraulic  jump.  The  Walnut  Creek  Flood-Control  Project  has  had 
much  attention  and  required  substantial  re-design.  As  with  most  urban  flood- 
control  channels,  increased  population  and  urbanization  generally  necessitates  a 
re-evaluation  of  the  design. 

This  research  can  be  broken  into  two  stages  as  illustrated  in  Figure  3.1.  The 
first  stage  deals  with  the  preliminary  issues  such  as  numerical  representation  of 
the  problem,  determination  of  the  design  criterion  or  objective  function,  and 
determination  of  the  design  variables  and  ranges  via  parameter  sensitivity  analysis. 
The  second  stage  deals  with  the  actual  optimization  of  the  problem.  There  must 
be  a  process  by  which  the  changes  to  the  design  variables  are  incorporated  into 
the  model.  Once  the  model  is  updated,  the  HIVEL2D  code  is  used  to  simulate  the 
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flow  for  the  new  design.  The  objective  function  is  evaluated  once  the  steady-state 
solution  is  reached,  and  the  design  space  derivatives  are  computed.  This  process  is 
repeated  for  a  prescribed  number  of  iterations  until  the  problem  has  converged  to 
an  answer.  Each  step  of  this  process  is  discussed  in  detail  in  the  following  sections. 

3.1  Numerical  Representation  of  Physical  Problem 

Representation  of  the  physical  problem  for  numerical  analysis  is  the  most 
important  part  of  numerical  modeling.  It  is  crucial  that  all  areas  vital  to  the 
behavior  of  the  system  are  well  represented,  thus  more  sophisticated  grids  are 
required. 

The  first  step  in  numerical  modeling  is  to  define  the  basic  geometry  of  the 
problem.  This  can  be  done  by  locating  important  points  in  the  design  such  as  the 
beginning  of  a  contraction  or  expansion.  Once  these  points  are  identified,  the  grid 
can  be  generated.  There  are  a  number  of  processors  for  generating  grids.  The  one 
used  for  this  work  is  the  Surface-Water  Modeling  System  or  SMS.  SMS  is  a  pre- 
and  post-processor  for  building  grids  and  viewing  solutions  [SMS  00]. 

Once  the  problem  geometry  is  represented  and  the  grid  is  generated,  the 
boundary  and  initial  conditions  can  be  specified.  At  this  point  the  problem  is  ready 
for  code  application  to  solve  for  the  flow  variables.  The  steady-state  solutions  can 
then  be  examined  with  a  post-processing  technique.  Using  the  results  from  the  flow 
simulation,  a  grid  convergence  test  is  conducted.  This  determines  the  appropriate 
grid  refinement  required  to  adequately  represent  the  problem. 
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Optimization 

Loop 


Figure  3.1:  Flow  Chart  of  Method  for  Optimizing  High-Velocity  Channel 
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3.2  Objective  Function 

A  variety  of  objectives  can  be  defined  for  high-velocity  channels.  The  objective 
functions  can  measure  variation  in  depth  to  produce  near  uniform  flow  depths, 
depth  raised  to  some  power,  and/or  energy  loss  through  a  transition.  The  functions 
can  be  functions  of  the  flow  variables,  the  grid,  and/or  the  design  variables.  A 
major  concern  for  engineers  designing  high-velocity  channels  is  water  overtopping 
the  walls  and  reaching  the  bridges.  To  avoid  this,  the  water  depth  and  cross-waves 
or  sloshing  should  be  minimized.  The  objective  function  used  by  Burg  is  given  by 

F0)  =  E  (jw(/3)  -  M/3))2  (3.1) 

where  h  is  the  depth  and  have  is  the  average  depth  over  the  entire  area  of 
interest  [Burg  99].  Equation  3.1  measures  the  non-uniformity  of  the  flow  depths 
by  measuring  the  variation  of  the  depths  from  the  average  depth  over  a  given 
area.  The  function  is  explicitly  dependant  only  on  Q  and  is  a  continuous  function; 
otherwise  the  design  space  derivatives  would  not  exist.  A  disadvantage  of  this 
function  is  its  inability  to  detect  subcritical  flow.  If  flow  becomes  subcritical,  the 
water  surface  is  smooth  thus  satisfying  the  equation,  but  the  depth  is  substantially 
increased. 

Another  possible  objective  function  is  given  by 

F0)  =  E  (M"  (3.2) 

Xk&^F 

where  n  equals  the  power.  Equation  3.2  is  the  sum  of  the  depth  raised  to  a  power. 
This  function  detects  increased  depths  such  as  would  occur  for  subcritical  flows. 
This  equation  is  not  as  sensitive  to  roughness  in  the  water  surface  as  Equation  3.1. 
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3.3  Parameter  Sensitivity  Analysis 

Several  model  parameters  exist  for  a  high-velocity  channel.  To  determine  the 
parameters  that  have  the  most  effect  on  the  model  globally  and  locally,  a  parameter 
sensitivity  analysis  is  performed.  Parameters  are  adjusted  independently,  and  the 
resulting  objective  function  values  are  recorded.  The  parameters  are  adjusted  to 
ensure  a  wide  range  of  possible  values  is  considered.  The  parameters  that  have 
the  greatest  effect  on  the  objective  function  are  determined.  Those  parameters 
that  have  little  or  no  effect  are  not  considered  in  the  optimization  process.  The 
parameters  that  point  to  an  obvious  solution  or  trend  through  the  sensitivity 
analysis  are  set  according  to  the  analysis  and  are  not  optimized.  The  resulting 
objective  function  values  are  also  used  to  determine  the  parameters  that  can  be 
optimized  locally  and  the  parameters  that  have  a  global  effect  on  the  model. 

3.4  Moving  the  Grid 

When  a  new  design  is  produced  by  the  optimization  code  the  new  grid  will 
need  to  be  developed  to  describe  the  new  design.  A  system  for  moving  the  mesh 
for  the  parameters  must  be  developed.  Moving  or  re-generating  the  grid  is  a  very 
important  part  of  the  process.  Most  of  an  engineer’s  time  is  devoted  to  generating 
an  accurate  representation  of  the  area  of  interest,  and  changes  to  the  mesh  must 
have  the  same  attention.  Three  ways  of  re-generating  a  grid  for  optimization 
are  integration  of  an  existing  grid  generation  code  within  the  optimization  code, 
develop  and  encode  a  set  of  rules  based  on  the  design  variables  and  grid  parameters, 
and  modify  an  existing  grid  based  on  the  design  variables. 

Most  grid  generation  tools  allow  users  to  save  steps  taken  to  generate  the 
grid.  These  steps  can  be  put  into  a  script  for  easy  re-generation  of  the  grid. 
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This  technique,  however,  requires  the  grid  to  be  re-generated  for  every  design 
variable.  For  small  problems,  re-generation  of  the  entire  grid  for  every  design 
variable  would  not  be  very  expensive.  The  re-generation  of  the  entire  grid  would 
be  quite  computationally  expensive  for  more  sophisticated  problems  of  larger  scale. 

Another  method  of  generating  the  new  grid  is  by  developing  and  encoding  a  set 
of  rules  based  on  the  design  variables  and  grid  parameters.  This  method  allows  for 
many  design  variables  and  grid  parameters  to  control  the  shape  of  the  boundary 
and  the  grid  spacing.  Once  the  boundary  is  determined,  a  grid  is  generated.  Burg 
[Burg  99]  uses  this  method  to  generate  structured,  multi-block  grids. 

The  final  method,  which  is  the  method  used  in  this  work,  modifies  an  existing 
grid  based  on  the  design  variables.  Due  to  the  fact  that  the  problems  used  in  this 
study  contain  several  hydraulic  structures,  a  structured  mesh  is  not  appropriate. 
The  mesh  used  is  an  unstructured,  two-dimensional  mesh.  Boundary  perturbations 
are  established  for  the  design  variables  using  the  optimization  procedure.  Once 
the  boundaries  have  been  moved,  a  call  is  made  to  a  Laplacian  solver  to  determine 
the  displacement  in  the  x  and  y  directions  of  the  interior  nodes  for  the  radius  of 
curvature  problem  or  in  the  z- value  for  the  super  elevation  problem  [Carey  99]. 
These  displacements  are  then  used  to  calculate  the  new  coordinates  for  the  nodes 
in  the  interior  of  the  domain.  The  equations  being  solved  are  given  by  Equation 
3.3.  This  technique  is  used  to  move  the  curve  for  the  various  radius  of  curvature 
values  and  is  illustrated  in  Figure  3.2.  This  technique  pushes  toward  a  uniform 
mesh  spacing  while  honoring  the  boundary  conditions  for  the  radius  of  curvature 
problem,  and  linearly  adjusts  the  z-values  for  the  super  elevation  problem. 


V2xd  =  0,  V2t/J  =  0,  V22d  =  0 


(3.3) 


36 


3.5  Flow  Simulation 

The  flow  analysis  is  performed  using  HIVEL2D.  HIVEL2D  is  a  finite  element 
model  that  solves  the  two-dimensional  shallow  water  equations.  The  results  are 
a  steady-state  representation  of  the  flow.  Details  of  HIVEL2D  are  given  in  the 
previous  chapter. 


3.6  Termination  Criteria 

The  optimization  terminates  when  the  termination  criteria  is  met.  For  this 
study,  the  optimization  criteria  is  given  by  Equation  3.4.  When  the  change  in  the 
objective  function  value  is  less  than  1CT7  the  code  terminates. 


F(pn)  -  F(pn-i)  <  10-7  (3.4) 

3.7  Function  Evaluation 

Once  the  steady-state  solution  is  obtained,  the  objective  function  is  evaluated. 
The  selected  objective  functions  are  expressed  in  the  non-linear  least  squares  form, 
thus  the  Gauss-Newton  optimization  algorithm  can  be  used.  That  is,  the  objective 
function  is  such  that  the  sum  of  the  squares  of  the  deviation  from  the  desired 
answer  is  minimized.  The  Gauss-Newton  method  takes  advantage  of  the  structure 
of  the  objective  function  to  approximate  the  Hessian.  Results  for  objective  function 
values  using  Equation  3.1  and  Equation  3.2  are  compared.  Generally,  the  equations 
give  similar  results.  Either  equation  can  be  used  and  similar  results  obtained.  The 
objective  function  used  in  this  study  is  given  in  Equation  3.1. 
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Figure  3.2:  Laplacian  Smoothing  Technique 
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3.8  Design  Space  Gradient  Calculation 

For  calculating  the  design  space  gradient,  the  fewer  steady-state  solution 
calculations  required  the  better.  This  is  due  to  the  fact  that  complex  problems,  as 
the  one  in  this  study,  contain  many  hydraulic  features  that  make  steady-state 
solution  calculations  computationally  expensive.  Sensitivity  analysis  does  not 
require  any  additional  steady-state  solutions  for  the  types  of  problems  addressed 
in  this  work  and  thus  is  quite  suitable  for  the  problems  discussed  in  this  work.  The 
discrete  sensitivity  analysis  requires  that  the  objective  function  be  evaluated  only 
at  steady  state;  therefore  the  gradient  estimation  routines  are  only  needed  once 
steady-state  flow  has  been  reached.  The  code  used  for  the  optimization  process 
gives  the  user  the  option  of  using  the  adjoint  variable  formulation  or  the  qnasi- 
analytic  method  to  calculate  the  derivatives.  The  adjoint  variable  formulation  is 
efficient  computationally  when  the  number  of  design  variables  is  larger  than  the 
number  of  objective  functions.  For  this  study,  the  goal  is  to  minimize  Equation 
3.1.  This  equation  is  in  the  general  non-linear  least  squares  form  illustrated  by 
Equation  3.5. 

N 

hs)  =  ;ea2  o?)  (3.5) 

k= 1 

fk(/3)  is  the  residual  function.  For  Equation  3.1,  fk{0)  —  have  ~  h.  In  this 
application,  there  is  one  residual  defined  at  each  of  the  41  nodes  illustrated  in 
Figure  3.3.  The  Gauss-Newton  method  is  applied  to  minimize  each  residual  with 
respect  to  the  selected  design  variable,  ffence,  for  the  problems  in  this  research 
there  are  41  constraint  equations  and  one  or  two  design  variables.  Consequently, 
the  quasi-analytic  formulation  would  be  the  most  efficient  method  of  calculating 
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the  necessary  design  space  gradients.  Using  this  formulation  enables  calculation  of 
the  design  space  gradients  by  solving  one  additional  linear  system  for  each  design 
variable  for  each  steady-state  HIVEL2D  simulation.  [Burg  99] 

3.9  Design  Variable  Update 

For  updating  the  design  variables,  Burg’s  code  gives  the  option  of  using  the 
method  steepest  descent  with  a  step-size,  the  method  of  steepest  descent  with 
a  linear  search,  the  Gauss-Newton  Algorithm,  or  the  BFGS  updating  method 
[Burg  99].  The  methods  used  in  this  research  are  the  method  of  steepest  descent 
with  a  step-size,  the  method  of  steepest  descent  with  a  linear  search,  and  the 
Gauss-Newton  Algorithm. 
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Figure  3.3:  Residuals  at  Each  Node  in  the  Objective  Space 


CHAPTER  IV 


RESULTS 

The  problem  at  Walnut  Creek  was  initially  addressed  using  a  physical 
model.  Engineers  using  the  existing  channel  and  the  engineering  design  manual 
recommendations  developed  an  initial  design.  The  design  was  then  used  to  build 
a  physical  model  of  the  site.  Once  the  physical  model  was  completed,  the  problem 
areas  and  potential  problem  areas  were  noted  and  modifications  were  designed. 
Adding  some  features,  such  as  the  divider  wall,  can  be  inexpensive.  To  test  various 
lengths  or  orientations  of  a  divider  wall,  sections  of  varying  lengths  are  constructed 
that  can  be  held  into  place  and  adjusted  for  orientation  by  engineers  until  the  flow 
conditions  improve.  Adjusting  curves  and  channel  widths  would  require  removing 
parts  of  the  model  and  rebuilding  them  with  the  new  design.  For  this  reason,  the 
physical  model  study  was  able  to  test  more  options  for  the  divider  wall  lengths 
and  orientation  than  for  the  other  design  parameters. 

The  first  step  to  apply  Burg’s  optimization  code  to  the  Walnut  Creek  study  is 
to  represent  the  site  numerically.  The  study  area  begins  at  the  junction  of  Walnut 
Creek  and  the  San  Ramon  Bypass  Channel  and  extends  upstream  in  the  San 
Ramon  Bypass  Channel  approximately  180  feet  as  indicated  in  Figure  4.1.  Through 
parameter  sensitivity  analysis  discussed  later  in  this  chapter,  the  super  elevation 
and  radius  of  curvature  in  the  San  Ramon  Bypass  Channel  (Figure  4.2)  are  selected 
as  the  system  design  variables  for  this  demonstration.  The  choice  of  an  objective 
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function  is  a  major  part  of  the  study.  For  an  optimization  problem,  the  objective 
function  must  contain  the  engineering  judgment.  All  of  the  necessary  requirements 
for  the  channel,  such  as  smooth  flow  and  minimal  depth,  must  be  included  in  the 
objective  function.  Using  the  objective  function,  a  parameter  sensitivity  analysis  is 
conducted  to  determine  the  design  parameters  that  have  an  affect  on  the  model.  It 
must  also  be  determined  whether  design  variables  have  global  or  local  effects.  The 
sensitivity  analysis  reduces  the  number  of  design  variables  to  only  those  to  which 
the  model  is  sensitive.  Once  the  design  parameters  are  determined,  Bing’s  code 
is  applied  to  find  an  improved  design.  The  design  parameters  for  this  study  are 
the  super  elevation  and  radius  of  curvature  of  the  curve  in  the  San  Ramon  Bypass 
Channel.  The  objective  is  to  minimize  overtopping  of  the  model  walls  and  water 
surface  roughness  just  upstream  of  the  confluence.  The  super  elevation,  denoted 
by  Ay  in  Figure  4.3,  and  radius  of  curvature,  denoted  by  R  in  Figure  4.4,  of  the 
curve  in  the  San  Ramon  Bypass  Channel  are  optimized  separately. 

4.1  Numerical  Representation  of  Physical  Problem 

The  first  step  in  generating  the  mesh  for  the  Walnut  Creek  study  is  to  identify 
important  points  in  AutoCAD.  Figure  4.5  illustrates  the  identification  of  important 
points  using  AutoCAD.  These  points  are  then  read  into  SMS  [SMS  00].  The  points 
are  assigned  a  bed  elevation  value  or  z- value  and  used  to  generate  a  number  of  nodes 
via  node  interpolation  using  linear  and  arc  interpolations.  Once  the  nodes  have 
been  created  the  triangulation  process  is  performed.  Once  the  mesh  is  generated, 
the  boundary  conditions  and  initial  conditions  for  the  Walnut  Creek  study  are 
applied. 


Figure  4.1:  The  Location  of  the  Study  Area. 
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Figure  4.3:  The  Super  Elevation  in  the  Curve. 


The  Walnut  Creek  problem  is  very  complex  and  includes  a  number  of  hydraulic 
phenomena.  For  this  reason,  the  steady-state  solution  is  more  difficult  to  obtain. 
In  order  to  prevent  introducing  artificial  hydraulic  jumps  that  are  hard  if  not 
impossible  for  the  model  to  push  out  of  the  system,  a  subcritical  tailwater  boundary 
was  set  to  maintain  a  reasonable  depth  throughout  the  model.  Once  the  jump  was 
pushed  through  the  model  to  the  outflow  boundary,  the  tailwater  boundary  was 
removed  and  the  outflow  was  set  to  supercritical. 

The  appropriate  grid  resolution  for  the  problem  is  determined  via  a  grid 
convergence  test.  This  is  accomplished  by  generating  an  initial  mesh  with  minimal 
resolution  and  increasing  the  resolution  until  the  changes  in  the  resulting  solutions 
are  minimal.  For  this  study,  the  areas  with  solutions  differing  by  more  than  10  3 
were  refined.  The  resulting  mesh  consists  of  4,592  elements  and  5,016  nodes.  This 
mesh  is  used  for  the  optimization  problems. 
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Figure  4.5:  AutoCAD  Depiction  of  Location  used  to  Generate  the  Mesh. 
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Once  the  mesh  is  created,  boundary  conditions  are  specihed.  The  Walnut 
Creek  problem  consists  of  three  boundary  conditions.  Inflow  boundaries  are  set 
for  Walnut  Creek  and  San  Ramon  Bypass  Channel.  The  inflow  for  Walnut  Creek 
is  supercritical  with  a  depth  of  11.83  feet  and  a  rate  of  6,800  cfs.  For  San  Ramon 
Bypass  Channel  the  inflow  is  supercritical  at  a  depth  of  15  feet  and  rate  of  15,200 
cfs.  These  values  represent  a  100-year  frequency  event.  An  outflow  boundary 
in  Walnut  Creek  is  set  to  supercritical  flow.  Figure  4.6  illustrates  the  mesh  and 
boundary  conditions.  The  Manning’s  n  value  for  the  entire  mesh  is  0.014. 

The  numerically  defined  problem  is  run  to  steady-state  using  HI\rEL2D.  The 
solution  is  compared  to  plots  given  in  the  physical  model  report  [Davis  87].  The 
comparisons  indicate  that  good  computational  fluid  dynamics  analysis  is  achieved. 
Figure  4.7  shows  a  qualitative  comparison  of  he  resulting  water  surface  elevations 
of  the  physical  model  and  the  numerical  model. 

4.2  Objective  Function 

The  Walnut  Creek  physical  model  study  considered  three  major  areas  when 
determining  an  appropriate  design.  The  areas  are  located  in  the  circular  curve 
section  of  Walnut  Creek  downstream  of  the  confluence,  in  the  section  of  Walnut 
Creek  located  just  upstream  of  the  confluence,  and  in  the  curved  section  of  the 
San  Ramon  Bypass  Channel  just  upstream  of  the  confluence.  However,  since  flow 
in  the  San  Ramon  Bypass  Channel  is  supercritical  everywhere,  it  is  only  affected 
by  changes  made  upstream  of  the  area.  For  this  reason,  we  will  only  consider  the 
first  two  areas  to  determine  what  has  an  affect  on  the  overall  model.  These  areas 
will  be  referred  to  as  area  1  and  area  2,  respectively,  and  are  shown  in  Figure 
4.8.  The  concern  in  the  areas  is  with  the  water  overtopping  the  walls.  This  can 
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Figure  4.7:  Water  Surface  Elevations  for  the  Physical  and  the  Numerical  Models 
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be  avoided  by  decreasing  the  overall  depth  and/or  by  decreasing  the  roughness 
of  the  water  surface.  Both  Equation  3.1  and  Equation  3.2  are  used  to  perform 
the  parameter  sensitivity  analysis.  The  parameter  sensitivity  analysis  discussed 
later  in  this  chapter  indicate  that  both  Equation  3.1  and  Equation  3.2  give  the 
same  results.  Further  tests  are  conducted  with  the  super  elevation  and  radius  of 
curvature  problems,  which  also  indicate  that  either  equation  could  be  used  and 
results  would  be  the  same.  For  this  study,  Equation  3.1  is  chosen  as  the  objective 
function  for  the  optimization  problems. 

4.3  Parameter  Sensitivity  Analysis 

In  a  practical  case,  such  as  this,  there  are  potentially  many  design  parameters. 
The  optimization  should  only  be  performed  on  the  design  parameters  that  have 
a  significant  impact  on  the  model.  A  parameter  sensitivity  analysis  is  used 
to  determine  the  design  variables  to  which  the  model  is  most  sensitive.  The 
primary  objectives  of  the  Walnut  Creek  physical  model  study  were  to  determine  the 
adequacy  of  the  San  Ramon  Bypass  Channel  and  the  Walnut  Creek- San  Ramon 
junction  and  to  develop  modifications  to  improve  the  adequacy.  The  major  areas 
considered  were  the  length  of  a  divider  wall  extension  located  at  the  junction,  the 
width  of  the  San  Ramon  Bypass  Channel,  and  the  super  elevation  and  radius  of 
curvature  of  the  curve  located  just  upstream  of  the  junction.  A  suite  of  test  cases 
is  run  to  assess  the  sensitivity  of  the  model  to  the  various  parameters.  The  list 
of  runs  made  to  test  the  model  sensitivity  are  given  in  Table  4.1.  Both  objective 
functions  were  evaluated  for  each  test  case  in  area  1  and  area  2. 

One  parameter  explored  is  the  super  elevation  in  the  curve  upstream  of  the 
confluence.  The  original  super  elevation  is  3.82  feet.  Three  other  super  elevation 
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Figure  4.8:  Objective  Space  for  Optimization. 
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Table  4.1:  Parameter  Sensitivity  Runs 


N  ame 

Width  of 

San  Ramon 

Channel  (ft) 

Super 

Elevation  (ft) 

Entrance 

Spiral 

Length  (ft) 

Exit 

Spiral 

Length  (ft) 

Radius  of 

Curve(ft) 

Divider 

Wall 

Length  (ft) 

Final  des 

23 

3.82 

200 

87.5 

285 

40 

SeO 

23 

0 

200 

87.5 

285 

40 

Se2 

23 

2 

200 

87.5 

285 

40 

Se6 

23 

6 

200 

87.5 

285 

40 

Width20 

20 

3.82 

200 

87.5 

285 

40 

Width25 

25 

3.82 

200 

87.5 

285 

40 

R350 

23 

3.82 

135.45 

49.72 

350 

40 

R400 

23 

3.82 

97.04 

9.997 

400 

40 

DwlO 

23 

3.82 

200 

87.5 

400 

10 

Dw30 

23 

3.82 

200 

87.5 

400 

30 

Dw64 

23 

3.82 

200 

87.5 

400 

63.8 

Crcl25 

23 

3.82 

0 

0 

125 

40 

Crc250 

23 

3.82 

0 

0 

250 

40 

Crc375 

23 

3.82 

0 

0 

375 

40 
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values  are  also  analyzed.  The  other  values  analyzed  are  no  super  elevation,  2  feet  of 
super  elevation,  and  6  feet  of  super  elevation.  Figure  4.9  shows  the  bed  elevations 
of  the  curve. 

The  resulting  objective  function  values  are  given  in  Table  4.2.  The  results 
indicate  that  the  model  is  not  very  sensitive  to  changes  to  the  super  elevation  in 
the  curve  in  area  1  and  area  2.  Figure  4.10  shows  the  water  surface  elevations  for 
the  original  design  and  the  three  test  cases  for  area  1  and  area  2.  Analyzing  the 
results  further  show  that  there  is  minimal  change  in  the  water  surface  elevations 
for  each  case.  The  resulting  depths  are  given  in  Figure  4.11.  The  water  surface 
elevation  values  are  differenced  for  each  of  the  three  cases.  The  resulting  differences 
for  the  water  surface  elevations  of  the  case  with  no  super  elevation  and  the  case 
with  a  super  elevation  value  of  2  feet  indicate  differences  of  approximately  0.15 
feet  occurring  in  a  few  locations  just  downstream  of  the  curve.  Examining  the 
differences  in  the  water  surface  elevations  for  the  case  with  no  super  elevation 
and  a  super  elevation  value  of  6  feet  indicate  differences  of  approximately  0.5 
feet  located  in  a  few  locations  just  downstream  of  the  bend.  Finally,  the  water 
surface  elevations  of  the  case  with  2  feet  of  super  elevation  and  6  feet  of  super 
elevation  are  compared.  Differences  are  visible  in  a  few  locations  just  downstream 
of  the  bend  with  values  of  approximately  0.25  feet.  Super  elevation  is  used  to 
suppress  the  disturbances  caused  by  curved  transitions,  which  arise  in  the  bend 
and  persist  downstream.  This  is  exactly  what  is  seen  from  the  cases  involving 
different  super  elevation  values.  The  affects  of  changing  the  super  elevation  appear 
in  the  downstream  end  of  the  bend  and  just  downstream  of  the  bend.  Area  1  and 
area  2  are  not  affected  by  these  changes.  This  indicates  that  the  super  elevation  in 
the  curve  can  be  optimized  without  considering  the  impact  on  area  1  and  area  2. 
In  order  to  accurately  optimize  for  super  elevation,  an  objective  function  should 


55 


Figure  4.9:  Bed  Elevations  for  No(left),  2’  (center),  and  6’  (right)  Super  Elevation. 
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be  located  in  the  area  impacted  by  changes  to  the  super  elevation.  In  this  case 
that  area  is  the  downstream  end  of  the  bend  and  just  downstream  of  the  bend. 


Table  4.2:  Super  Elevation  Parameter  Sensitivity  Objective  Function  Values 


Run  Name 

Objective  Function  1 

Objective  Function  2 

Areal 

Area2 

Areal 

Area2 

SeO 

1.140224 

1.519831 

12589.39 

5344.084 

Se2 

1.141895 

1.528477 

12589.71 

5349.77 

Se6 

1.14185 

1.637631 

12597.34 

5357.59 

The  radius  of  curvature  in  this  curve  is  also  adjusted.  The  curve  is  restricted 
to  a  section  of  the  channel  that  would  not  require  additional  changes  to  the  design. 
That  is,  given  two  fixed  channels  design  a  curved  transition  including  entrance  and 
exit  spirals  to  join  the  two.  This  restraint  limits  the  radius  of  curvature  to  values 
between  285  feet  and  410  feet.  The  initial  radius  of  curvature  is  285  feet.  The 
lengths  of  the  entrance  and  exit  spirals  are  200  feet  and  87.5  feet,  respectively. 
The  different  values  of  the  radius  of  curvature  tested  were  a  radius  of  curvature 
of  350  feet  with  corresponding  entrance  and  exit  spiral  lengths  of  135.45  feet  and 
49.72  feet,  respectively,  and  a  radius  of  curvature  of  400  feet  with  corresponding 
entrance  and  exit  spiral  lengths  of  97.04  feet  and  9.997  feet,  respectively.  The 
layout  of  the  channels  including  bed  elevations  are  depicted  in  Figure  4.12. 

The  steady-state  solutions  for  the  different  scenarios  indicate  that  area  1  and 
area  2  are  not  affected  by  changes  in  the  radius  of  curvature.  This  is  shown  by  the 
objective  functions  given  in  Table  4.3.  The  water  surface  profiles  in  area  1  and  area 
2  of  the  radius  of  curvature  test  cases,  illustrated  in  Figure  4.13,  shows  that  the 
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Figure  4.10:  Depth  Profile  in  Areal  and  Area2  for  Super  Elevation  Test  Cases. 
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Figure  4.11:  Resulting  Depths  for  the  Super  Elevation  Test  Cases. 
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changes  in  the  surface  water  elevations  is  negligible  for  each  change  in  the  radius 
of  curvature.  The  impact  of  the  changes  to  the  curved  transition  appear  in  the 
curve  and  just  downstream  of  the  curve.  This  can  be  optimized  for  local  concerns 
without  regard  to  the  area  1  and  area  2  impacts.  For  this  study,  the  objective 
function  is  located  within  the  bend  and  just  downstream  of  the  bend  to  optimize 


the  radius  of  curvature. 


^  «»  ^  ^ 


Figure  4.12:  Bed  Elevations  of  the  Original  Curve,  350’  Radius,  and  400’  Radius. 
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Figure  4.13:  Depth  Profiles  in  Areal  and  Area2  for  Radius  of  Curvature. 
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Table  4.3:  Radius  of  Curvature  Parameter  Sensitivity  Objective  Function  Values 


Run  Name 

Objective  Function  1 

Objective  Function  2 

Areal 

Area2 

Areal 

Area2 

Finaldes 

1.145559 

1.576477 

12586.15 

5357.123 

R350 

1.136146 

1.604369 

12585.86 

5354.309 

R400 

1.151916 

1.835092 

12583.6 

5374.444 

A  third  area  examined  deals  with  the  width  of  the  San  Ramon  Channel.  This 
channel  has  an  initial  width  of  32  feet  upstream  and  contracts  to  23  feet  upstream 
from  the  junction  with  Walnut  Creek.  Widths  of  20  feet  and  25  feet  are  tested.  The 
widths  are  adjusting  by  moving  the  wall  farthest  from  Walnut  Creek  only.  This  is 
done  to  avoid  introducing  any  adverse  affects  from  changing  the  divider  wall.  The 
transition  from  32  feet  to  the  new  width  became  unsymmetrical  as  a  result,  but 
the  affects  were  negligible.  Figure  4.14  shows  the  mesh  outline  and  bed  elevations 
for  the  two  test  cases.  The  resulting  objective  function  values  are  given  in  Table 
4.4.  Figure  4.15  shows  the  resulting  water  surface  elevations  for  channel  width  of 
20  feet,  25  feet,  and  the  original  width  of  23  feet.  The  resulting  depths  are  shown 
in  Figure  4.16.  The  objective  functions  indicate  that  the  narrower  the  channel  the 
better  the  design.  However,  with  a  narrower  channel  come  higher  depths.  The 
wall  heights  in  San  Ramon  are  15  to  22  feet.  The  depths  in  San  Ramon  Bypass 
Channel  for  the  20-foot  wide  problem  are  too  near  the  top  of  the  walls.  Without 
knowing  the  cost  of  moving  sections  of  the  wall,  the  decision  is  made  to  use  the 
results  from  the  physical  model  study,  that  is,  a  channel  width  of  23  feet. 
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Figure  4.14:  Bed  Elevations  for  Channel  Widths  of  20’  (left)  and  25’  (right). 
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Figure  4.15:  Water  Depth  Profiles  in  Areal  and  Area2  for  Width  Test  Cases. 
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Figure  4.16:  Depths  for  Channel  Widths  of  20’  (left),  23’  (center),  and  25’  (right). 
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Table  4.4:  Channel  Width  Parameter  Sensitivity  Objective  Function  Values 


Run  Name 

Objective  Function  1 

Objective  Function  2 

Areal 

Area2 

Areal 

Area2 

Width20 

1.335859 

0.9924769 

12203.65 

4991.545 

Width23 

1.145559 

1.576477 

12586.15 

5357.123 

Width25 

1.284566 

3.663777 

12946.4 

5610.11 

The  final  area  explored  is  the  length  of  the  divider  wall  extension  located  at  the 
confluence.  Test  with  lengths  of  the  divider  wall  at  10  feet,  30  feet,  40  feet,  and  63.8 
feet  were  performed.  Table  4.5  shows  that  the  objective  functions  are  sensitive  to 
the  changes  in  the  divider  wall  extension  lengths.  The  values  of  objective  function  1 
in  area  1  indicate  that  the  longer  divider  wall  lengths  are  better.  However,  objective 
function  2  in  area  2  indicates  just  the  opposite.  Figure  4.17  shows  the  resulting 
water  surface  profiles  for  the  divider  wall  test  cases.  Examination  of  the  resulting 
depths  reveals  the  reasoning.  It  is  possible  to  generate  a  smooth  water  surface 
elevation  while  increasing  the  depth.  Subcritical  flow  has  an  increased  depth  but 
a  smoother  water  surface.  While  using  an  objective  function  that  considers  the 
roughness  of  the  water  surface,  it  is  possible  to  overlook  increased  depth.  The 
results  show  a  hydraulic  jump  located  in  area  2.  As  the  divider  wall  is  lengthened, 
the  jump  spreads  further  upstream.  The  jump  eventually  spreads  throughout  area 
2  resulting  in  a  smoother  water  surface.  Overtopping  occurs,  however,  because 
the  depth  is  increased.  The  length  of  the  divider  wall  can  be  determined  using 
sensitivity  analysis.  For  this  reason  the  divider  wall  length  determined  by  the 
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physical  model  study  is  used  for  this  study,  and  the  length  of  the  divider  wall  is 
not  optimized. 


Table  4.5:  Divider  Wall  Length  Parameter  Sensitivity  Objective  Fimction  Values 


Run  Name 

Objective  Function  1 

Objective  Function  2 

Areal 

Area2 

Areal 

Area2 

DW10 

3.27303652 

429.7831586 

19958.786 

7223.307 

DW30 

3.449422665 

284.275526 

19652.89796 

12108.97285 

DW40 

3.212145535 

178.1307028 

19694.6363 

12442.91466 

DW64 

6.664623238 

4.86675199 

19289.05433 

14005.4 

The  parameter  sensitivity  analysis  indicate  that  the  parameters  that  have  the 
most  impact  on  Areal  and  Area2  are  the  width  of  San  Ramon  Bypass  Channel 
and  the  length  of  the  divider  wall  extension.  The  super  elevation  and  the  radius 
of  curvature  of  the  curve  section  in  San  Ramon  Bypass  Channel  have  a  minimal 
impact  on  the  model  globally.  However,  the  results  from  the  sensitivity  analysis 
for  the  width  of  San  Ramon  Bypass  Channel  indicate  that  the  channel  should 
be  as  narrow  as  possible  without  going  subcritical  to  achieve  a  better  design. 
Narrower  channels,  however,  lead  to  higher  depths  in  San  Ramon  Bypass  Channel 
and  potential  overtopping.  The  depths  for  the  channel  width  of  20  feet  are  too  close 
to  the  top  of  the  walls.  To  ensure  that  an  appropriate  factor  of  safety  is  met,  the 
channel  width  of  23  feet  is  used.  Likewise,  the  divider  wall  sensitivity  study  also 
points  to  an  answer.  The  longer  the  divider  wall,  the  better  the  flow  downstream 
from  the  confluence.  However,  if  the  wall  is  too  long,  a  hydraulic  jump  is  pushed 
upstream  which  eventually  overtops  the  channel  wall.  The  suggested  design  is 
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Figure  4.17:  Water  Depth  Profiles  in  Areal  and  Area2  for  Divider  Wall  Test  Cases. 
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to  lengthen  the  divider  wall  until  the  hydraulic  jump  is  within  a  given  distance 
from  the  shorter  channel  wall.  Though  the  results  for  the  sensitivity  analysis  of 
the  length  of  the  divider  wall  has  an  effect  globally  on  the  model,  without  more 
information  about  constraints  in  the  area  like  factor  of  safety  with  regards  to  the 
wall  heights  it  is  not  possible  to  determine  a  more  adequate  design  than  that 
determined  by  the  physical  model  study. 

The  model  is  locally  sensitive  to  the  super  elevation  and  the  radius  of  curvature 
of  the  curved  section  in  San  Ramon  Bypass  Channel.  That  is,  changes  to  the  super 
elevation  and  the  radius  of  curvature  create  great  disturbances  in  the  curved  section 
in  San  Ramon  Bypass  Channel.  These  parameters  can  be  optimized  locally.  That 
is  the  design  considerations  for  these  two  parameters  are  in  the  San  Ramon  Bypass 
Channel  only.  This  work  uses  optimization  techniques  to  determine  a  better  design 
for  the  super  elevation  and  radius  of  curvature  of  the  curve  in  the  San  Ramon 
Bypass  Channel.  The  results  are  discussed  in  the  following  chapter. 

4.4  Moving  the  Grid 

A  Laplacian  smoothing  technique,  illustrated  in  the  previous  chapter,  is  used 
to  smooth  the  grid  once  the  boundaries  have  been  moved.  The  difficulty  is  in 
locating  the  boundaries  for  the  given  changes  to  the  design  parameter.  For  the 
super  elevation  problem,  the  z-values  are  adjusted  by  calculating  a  percent  change 
for  a  given  length.  Once  the  boundaries  are  set,  the  remaining  z-values  are  adjusted 
using  the  Laplacian  technique.  The  entrance  and  exit  spirals  complicate  adjusting 
the  boundaries  for  the  radius  of  curvature  problem.  For  the  radius  of  curvature 
problem,  the  tangent  lengths  remain  constant  to  insure  that  the  mesh  is  only  moved 
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in  the  area  of  the  current  curve.  With  a  given  value  for  the  radius  of  curvature,  the 
lengths  of  the  entrance  and  exit  spirals  are  obtained  using  the  following  equations. 


T\  =  X\  —  R  *  sinAi  + 


Y2  +  R*  cosA2 


(Y\  +  R*  cos  Ai)  *  cos  I 
sinl 


(4.1) 


T2 


X2  —  R*  sinA2  + 


Y\  +  R  *  cosAi  —  (Y2  +  R  *  cosA2 )  *  cosl 
sinl 


(4.2) 


where  Tj  and  T2  are  the  tangent  lengths  for  the  entrance  and  exit  spirals, 
respectively,  R  is  the  radius  of  curvature,  X\  and  Y\  are  the  coordinates  of  the 
point  of  change  from  the  entrance  spiral  to  the  curve,  X2  and  Y2  are  the  coordinates 
of  the  point  of  change  from  the  curve  to  the  exit  spiral,  Ai  is  the  central  angle  of 
the  entrance  spiral,  and  A2  is  the  central  angle  of  the  exit  spiral  [Rubey  38].  The 
equations  are  solved  using  a  Newton  iterative  solver.  A  layout  of  the  curve  is  given 
in  Figure  4.18. 


4.5  Circular  Curve  Problem 

To  test  the  ability  of  the  model  to  produce  reasonable  answers,  a  test  case  is 
conducted.  The  radius  of  a  curve  with  no  spirals  or  super  elevation  is  expected  to 
increase  to  generate  smoother  flow  downstream.  A  model  representing  the  Walnut 
Creek  area  is  generated  with  a  circular  curve  with  no  entrance  or  exit  spirals 
and  no  super  elevation  in  the  San  Ramon  Bypass  Channel  just  upstream  of  the 
confluence.  The  initial  design  is  illustrated  in  Figure  4.19.  For  this  design  problem, 
the  radius  of  the  circular  curve  is  the  only  design  variable.  The  objective  function 
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is  a  measure  of  the  roughness  or  non- uniformity  of  the  water  surface  in  the  curve. 
Figure  4.20  shows  the  location  of  the  objective  function  and  the  design  space.  The 
design  space  consists  of  a  circular  curve  with  two  straight  sections.  The  expected 
optimum  would  be  to  push  the  curve  to  the  design  space  limits  and  eliminate 
the  straight  sections.  This  would  create  a  smoother  transition,  thus  reducing  the 
movement  of  the  water. 

The  radius  of  the  curve  is  set  at  125  feet  initially  and  allowed  to  move  between 
0  feet  and  400  feet.  The  points  of  tangency  with  the  model  are  allowed  to  move 
but  are  constrained.  Figure  4.22  shows  the  points  of  tangency  for  the  maximum 
radius  of  curvature  value.  The  design  process  terminates  when  the  design  variable 
(radius)  becomes  larger  than  the  constraint  of  400  feet.  There  is  an  upper  bound 
because  the  end  points  are  fixed.  The  optimization  parameters  consist  of  using 
the  quasi-analytic  method  to  compute  the  derivatives  and  the  method  of  steepest 
descent  with  a  linear  search  to  update  the  design  variable.  Figure  4.21  shows  the 
objective  function  values  for  the  corresponding  radius  values.  The  results  show  that 
the  larger  the  radius  the  better  the  design.  This  indicates  the  need  for  entrance 
and  exit  spirals  and  super  elevation. 

4.6  Super  Elevation 

When  changing  the  direction  of  flow  via  curved  channels,  the  centripetal  force 
of  the  fluid  causes  a  rise  in  the  water  surface  on  the  outside  wall  and  a  depression 
in  the  water  surface  along  the  inside  wall.  These  curves  also  cause  disturbances 
in  the  flow  that  can  persist  downstream.  To  balance  the  force  of  the  fluid,  the 
channel  bed  can  be  sloped  across  the  channel.  This  transverse  bed  slope  is  known 
as  super  elevation. 
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Figure  4.19:  Initial  Mesh  for  the  Circular  Curve  Optimization  Problem. 
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Figure  4.21:  Objective  Function  Values  versus  Radius  of  Curvature. 
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Figure  4.22:  Layout  of  Curve  Including  Bounds  for  Points  of  Tangency. 
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The  curved  section  of  the  San  Ramon  Bypass  Channel,  shown  in  Figure  4.23, 
changes  the  flow  direction  by  43.6  degrees.  The  flow  is  supercritical  and  at  a  rate 
of  15,200  cfs.  The  channel  is  23  feet  wide  in  the  curved  section.  The  suggested 
range  of  super  elevation  values  for  these  flow  and  geometry  parameters  given  by 
Equation  4.3  is  between  2.137906202  feet  and  4.2758124  feet  [EM  1110-2-1601]. 
This  is  based  on  the  assumption  that  the  velocity  distribution  is  uniform  across  the 
channel  and  that  the  radius  of  curvature  is  285  feet  for  the  entire  curve.  However, 
the  velocity  distribution  is  not  likely  to  be  uniform  across  the  channel  and  the 
radius  of  curvature  varies  from  273.5  feet  and  296.5  feet  across  the  channel.  The 
super  elevation  determined  by  the  physical  model  study  was  3.82  feet. 


A  y  =  C 


V2W 

gr 


(4.3) 


where 

Ay  =  super  elevation 

C  =  coefficient  between  0.5  and  1.0 

V  =  mean  channel  velocity 

W  =  channel  width 

g  =  acceleration  of  gravity 

r  =  radius  of  curvature  to  center  line 

The  design  variable  is  the  super  elevation  in  the  curve  in  San  Ramon 
Bypass  Channel  (Figure  4.24).  The  problem  is  optimized  using  Equation  3.1 
as  the  objective  function.  The  objective  space  is  located  along  the  walls  in  the 
downstream  section  of  the  curve  and  just  downstream  of  the  curve  as  illustrated  in 
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Figure  4.23:  Initial  Grid  for  the  Super  Elevation  Problem. 
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Figure  4.25.  The  initial  design  has  a  super  elevation  value  of  2.14  feet.  The  design 
iterations  are  plotted  against  the  design  variable  in  Figure  4.26.  The  improved 
solution  is  determined  within  nine  design  iterations.  The  resulting  value  for  super 
elevation  is  3.35  feet  and  is  illustrated  in  Figure  4.27.  Comparing  this  value  to 
the  value  determined  by  the  physical  model,  which  is  3.82  feet,  and  the  range 
determined  by  Equation  4.3,  which  is  2.138  feet  to  4.28  feet,  shows  that  the  result 
from  the  optimization  code  is  reasonable.  Furthermore,  comparing  the  depth 
profiles  for  the  initial  design  and  new  design  (Figure  4.28)  shows  that  the  new 
design  is  an  improved  design.  The  maximum  water  depth  is  reduced  by  0.3  feet.  In 
some  locations  the  depth  is  reduced  by  more  than  0.5  feet.  The  objective  function 
value  is  reduced  in  the  improved  design  by  30%.  Figure  4.29  shows  the  depth 
profile  of  the  physical  model  design  and  the  improved  design.  The  profiles  are 
very  similar  and  indicate  the  ability  of  the  optimization  code  to  produce  improved 
channel  design.  A  plot  of  the  design  variable  versus  the  objective  function  is  given 
in  Figure  4.30. 


4.7  Radius  of  Curvature 

Curved  channels  are  used  to  change  the  direction  of  flow.  To  transition  the 
flow  more  smoothly,  entrance  and  exit  spirals  are  used  to  gradually  change  the 
direction  of  flow.  A  poorly  designed  curve  can  cause  disturbances  in  the  flow 
that  can  persist  downstream.  The  proper  combination  of  entrance  and  exit  spiral 
lengths  and  radius  of  curvatures  for  the  circular  curve  are  required  to  insure  the 
least  amount  of  disturbances  possible. 

The  curve  section  of  the  San  Ramon  Bypass  Channel,  shown  in  Figure  4.39, 
changes  the  flow  direction  by  43.6  degrees.  The  flow  is  supercritical  and  at  a 


Figure  4.25:  The  Objective  Area  for  the  Super  Elevation  Problem. 
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Figure  4.27:  Layout  and  Bed  Elevations  of  the  Improved  Design. 
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—  Initial  Design  Outer  Wall 

—  Initial  Design  Inner  Wall 
Improved  Design  Outer  Wall 
Improved  Design  InnerWall 


Station  Location  (ft) 


Figure  4.28:  Depth  Profiles  for  Initial  Design  and  Improved  Design. 
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—  Physical  Model  Design  Outer  Wall 

—  Physical  Model  Design  InnerWall 
Improved  Design  Outer  Wall 
Improved  Design  InnerWall 
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Figure  4.29:  Depth  Profiles  for  Physical  Model  Design  and  Improved  Design. 
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Super  Elevation  Value 


Figure  4.30:  Design  Variable  versus  Objective  Function 
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rate  of  15,200  cfe.  The  channel  is  23  feet  wide  in  the  curved  section.  The  design 
determined  by  the  physical  model  study  gives  a  radius  of  curvature  of  285  feet, 
entrance  spiral  length  of  200  feet  with  a  tangent  length  of  209.02  feet,  and  an  exit 
spiral  length  of  87.5  feet  with  a  tangent  length  of  164.99  feet  [Davis  87]. 

The  design  variable  is  the  radius  of  curvature  in  the  curve  in  San  Ramon 
Bypass  Channel  (Figure  4.32).  The  problem  is  optimized  using  Equation  3.1 
as  the  objective  function.  The  objective  space  is  located  along  the  walls  in  the 
downstream  section  of  the  curve  and  just  downstream  of  the  curve  as  illustrated 
in  Figure  4.33.  The  initial  design  has  a  radius  of  curvature  of  325  feet.  The  design 
iterations  are  plotted  against  the  design  variable  in  Figure  4.34.  The  improved 
solution  is  determined  within  30  design  iterations.  The  resulting  value  for  radius 
of  curvature  is  302  feet  and  is  illustrated  in  Figure  4.35.  Comparing  this  value 
to  the  value  determined  by  the  physical  model,  which  is  285  feet,  shows  that 
the  result  from  the  optimization  code  is  reasonable.  Furthermore,  comparing  the 
depth  profiles  for  the  initial  design  and  new  design  (Figure  4.36)  shows  that  the 
new  design  is  an  improved  design.  The  maximum  water  depth  is  reduced  by  0.4 
feet.  In  some  locations  the  depth  is  reduced  by  more  than  0.5  feet.  The  objective 
function  value  is  reduced  by  17%.  Figure  4.37  shows  the  depth  profile  of  the 
physical  model  design  and  the  improved  design.  The  profiles  are  very  similar  and 
indicate  the  ability  of  the  optimization  code  to  produce  improved  channel  design. 
A  plot  of  the  design  variable  versus  the  objective  function  is  given  in  Figure  4.38. 

4.8  Radius  of  Curvature  and  Super  Elevation 

There  are  several  design  parameters  for  a  high-velocity  channel.  The  ideal 
design  improvement  code  would  consider  all  of  these  parameters  when  determining 
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Figure  4.35:  Layout  of  the  Improved  Design. 
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Figure  4.36:  Depth  Profiles  for  Initial  Design  and  Improved  Design. 


93 


Figure  4.37:  Depth  Profiles  for  Physical  Model  Design  and  Improved  Design. 
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a  design.  To  determine  the  ability  of  the  code  to  handle  multiple  design  variables, 
a  test  was  conducted  to  assess  the  codes  ability  to  optimize  on  the  radius  of 
curvature  and  super  elevation  simultaneously.  Figure  4.40  shows  the  design  space 
for  the  two  design  variable  problem.  The  problem  is  optimized  using  Equation  3.1 
as  the  objective  function.  The  objective  space  is  located  along  the  walls  in  the 
downstream  section  of  the  curve  and  just  downstream  of  the  curve  as  illustrated 
in  Figure  4.41.  The  design  iterations  are  plotted  against  the  design  variable  in 
Figure  4.42.  The  improved  solution  is  determined  within  20  design  iterations.  The 
resulting  value  for  radius  of  curvature  is  301.9  feet  and  is  illustrated  in  Figure  4.43. 
Comparing  this  value  to  the  value  determined  by  the  physical  model,  which  is  285 
feet,  shows  that  the  result  from  the  optimization  code  is  reasonable.  Furthermore, 
comparing  the  depth  profiles  for  the  initial  design  and  new  design  (Figure  4.44) 
shows  that  the  new  design  is  an  improved  design.  The  maximum  water  depth  is 
reduced  by  0.5  feet.  In  some  locations  the  depth  is  reduced  by  more  than  0.85 
feet.  The  objective  function  value  is  reduced  by  40%.  Figure  4.45  shows  the  depth 
profile  of  the  physical  model  design  and  the  improved  design.  The  profiles  are 
very  similar  and  indicate  the  ability  of  the  optimization  code  to  produce  improved 
channel  design.  A  plot  of  the  design  variable  versus  the  objective  function  is  given 
in  Figure  4.46,  Figure  4.47,  and  Figure  4.48. 
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Figure  4.43:  Layout  of  the  Improved  Design. 
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Figure  4.44:  Depth  Profiles  for  Initial  Design  and  Improved  Design. 
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Figure  4.45:  Depth  Profiles  for  Physical  Model  Design  and  Improved  Design. 
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Figure  4.46:  Design  Variable  versus  Objective  Function. 
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Design  Variable  versus  Objective  Function. 
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Figure  4.48:  Design  Variable  versus  Objective  Function. 


CHAPTER  V 


SUMMARY,  CONCLUSIONS,  AND  RECOMMENDATIONS  FOR  FUTURE 

WORK 

5.1  Summary 

The  capability  to  use  optimization  techniques  to  design  “real-world”  high- 
velocity  channels  has  been  tested,  and  the  methodology  for  such  an  application 
has  been  developed.  The  results  show  that  optimization  techniques  along  with  an 
automated  process  to  re-generate  the  grid  are  capable  of  determining  an  improved 
design  for  high-velocity  channels. 

Several  preliminary  steps  are  necessary  before  the  optimization  techniques  are 
applied.  Developing  a  numerical  representation  of  the  physical  problem  is  the 
first  preliminary  step.  After  the  physical  problem  is  represented  numerically, 
an  objective  function  and  the  objective  space  are  determined.  This  is  a  crucial 
part  of  the  optimization  process,  since  the  objective  function  must  encompass  all 
of  the  engineering  judgment.  The  objective  function  in  this  work  measures  the 
non-uniformity  of  the  flow  depths.  High-velocity  channels  consist  of  many  design 
parameters.  Some  of  these  parameters  have  more  of  an  impact  on  the  channel  than 
others.  The  parameters  to  which  the  model  is  most  sensitive  can  be  determined  via 
parameter  sensitivity  analysis.  This  allows  for  limitation  of  design  parameters  to 
those  with  the  greatest  impact  on  the  channel  and  indicates  whether  sensitivities 
to  the  parameters  are  global  or  local. 
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The  application  of  the  optimization  techniques  to  a  high-velocity  channel 
requires  a  mechanism  to  reflect  changes  to  the  design  variables  in  the  channel. 
This  is  accomplished  through  an  automated  grid  generation  tool.  Once  the 
channel  design  is  updated  with  respect  to  the  design  variables,  HIVEL2D  is  used 
to  simulate  the  steady-state  flow  conditions.  These  steady-state  conditions  are 
used  to  determine  the  objective  function  value.  The  objective  function  value 
is  indicative  of  the  adequacy  of  the  design.  The  quasi- analytic  formulation  is 
used  to  compute  the  design  space  gradients,  and  the  design  variables  are  updated 
using  the  Gauss-Newton  method.  The  optimization  process  is  repeated  until  the 
model  is  converged.  The  convergence  criterion  is  based  on  the  change  in  the 
objective  function  for  the  current  design  iteration  and  the  previous  design  iteration. 
Convergence  of  the  optimization  code  yields  an  improved  design. 

The  super  elevation  and  radius  of  curvature  in  the  bend  just  upstream  of  the 
Walnut  Creek-San  Ramon  Bypass  Channel  junction  are  determined  using  this 
optimization  technique.  Each  design  parameter  has  a  local  effect.  That  is,  changes 
to  the  super  elevation  and  radius  of  curvature  in  the  bend  effect  the  model  in 
the  bend  and  just  downstream  of  the  bend  and  have  minimal  impact  in  other 
areas  of  the  channel.  The  optimization  process  is  applied  to  each  design  variable 
separately,  and  then  it  is  applied  to  both  design  parameters  via  two  design  variable 
optimization. 


5.2  Conclusions 


The  application  of  the  optimization  process  to  the  design  variables  of  super 
elevation  and  radius  of  curvature  in  the  bend  of  San  Ramon  Bypass  Channel 
yields  an  improved  design.  The  objective  function  for  the  super  elevation  problem 
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is  reduced  by  30%.  A  reduction  of  the  objective  function  of  17%  is  achieved  in 
the  radius  of  curvature  design  problem.  By  applying  the  optimization  technique 
to  both  super  elevation  and  radius  of  curvature  through  a  multiple  design  variable 
process,  it  is  possible  to  reduce  the  objective  function  value  by  over  40%.  As 
indicated  by  the  objective  function,  each  case  provides  an  improved  channel  design. 

The  results  show  that  optimization  techniques  can  be  applied  to  realistic  high- 
velocity  channels  to  produce  an  improved  design.  The  steps  by  which  such  an 
application  is  made  are  developed  and  outlined  in  this  work.  The  approach 
has  much  utility  in  high-velocity  channel  design.  Using  this  technique  will  allow 
screening  of  designs  for  numerous  design  parameters.  With  a  more  sophisticated 
automation  process,  it  would  be  possible  to  design  an  entire  channel  to  meet 
certain  constraints  and  produce  a  prescribed  flow  condition.  Although  this  research 
addressed  a  two-dimensional  problem  that  solved  the  shallow  water  equations,  the 
optimization  techniques  are  not  limited  to  such  problems.  They  can  also  be  applied 
to  various  types  of  flow  problems  that  solve  the  Navier-Stokes  equations  and  to 
groundwater  flow  problems,  as  well  as  problems  in  other  engineering  disciplines. 

5.3  Recommendations  for  Future  Work 

There  are  a  few  areas  that  can  be  investigated  to  expand  this  technique. 
Development  of  an  improved  automation  process  is  important  for  expanding  the 
process  to  different  problems  and  disciplines.  The  current  automation  process 
is  designed  for  the  problem  investigated  in  this  work.  Some  effort  is  required 
to  apply  this  technique  to  super  elevation  and  radius  of  curvature  in  another 
problem.  Applying  this  technique  to  other  design  parameters  would  require  a 
re-evaluation  of  the  automation  process.  The  automation  process  that  would  be 
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most  valuable  would  be  one  that  has  the  capability  of  generating  the  entire  mesh 
given  prescribed  parameters  such  as  fixed  boundaries,  moveable  boundaries,  and 
boundary  conditions.  This  could  be  enhanced  via  a  graphical  user  interface,  which 
allowed  the  problem  to  be  defined,  launch  the  optimization  technique,  and  provide 
real-time  updates  of  the  design  variables.  Such  a  process  would  make  it  easier  to 
apply  the  optimization  techniques  to  multiple  design  variables. 

The  problem  examined  in  this  work  dealt  with  the  two-dimensional  shallow 
water  equations.  Additional  work  should  be  done  to  apply  these  optimization 
strategies  to  three-dimensional  problems,  thus  allowing  investigation  of  a  wider 
range  of  channels.  An  important  issue  to  remember  when  applying  the 
optimization  techniques  to  other  problems  is  that  the  objective  function  is  site- 
specific,  and  the  design  criteria  are  different  for  every  problem.  Because  of  these 
differences,  the  objective  function  should  be  investigated  for  every  application. 
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Engineers  are  sometimes  required  to  select  the  “best”  design  from  a  set  of 
feasible  designs.  This  is  called  optimal  design  or  optimization.  Optimization 
principles  are  of  increasing  importance  in  modern  design.  Even  in  standard 
applications,  much  is  to  be  gained  by  using  an  algorithm  to  generate  initial  meshes 
and  refine  to  a  “near-optimal”  mesh  automatically. 

For  this  work  optimization  is  generally  defined  as  follows.  Given  a  function 
of  one  or  more  independent  variables,  End  the  value  of  those  variables  where  F  is 
a  minimum.  The  function  F  should  give  an  objective  measure  of  the  quality  of 
the  system  state.  Minimizing  F  should  be  accomplished  using  as  few  calculations 
of  F  as  possible  due  to  the  propensity  of  the  calculations  to  be  computationally 
expensive.  Values  of  the  variables  are  calculated  using  iterative  processes  that  start 
at  some  point  and  move  stepwise  to  points  for  which  F  is  smaller.  There  are  many 
different  means  of  calculating  the  next  guess  for  the  variables.  Two  major  methods 
are  heuristic  methods  and  gradient-based  methods.  The  influence  coefficients  are 
calculated  via  sensitivity  analysis. 

Heuristic  Method 

Heuristic  search  techniques,  such  as  genetic  algorithms,  simulated  annealing, 
and  neural  networks  achieve  nearly  optimal  solutions  at  a  reasonable  computational 
cost,  without  guaranteeing  a  globally  optimal  solution  will  be  found.  These 
techniques  are  stochastic  search  procedures  that  use  probabilistic  rather  than 
deterministic  search  rules.  Genetic  algorithms  are  search  algorithms  based  on 
the  concept  of  natural  selection  and  natural  genetics  [Holland  75].  The  objective 
function  magnitude,  instead  of  derivative  information,  is  used  directly  in  the  search, 
therefore  allowing  these  techniques  to  be  applied  to  nonconvex,  highly  nonlinear 


117 


and  complex  problems  [Goldberg  89].  Techniques  of  this  type  can  locate  near- 
optimal  solutions  of  complex  problems  with  very  large  nonlinear  search  spaces. 
Discontinuous  fitness  functions  and  discrete  input  variables  may  be  used.  Genetic 
algorithms  do  not  fall  into  local  minima  easily,  and  the  configuration  decisions 
proceed  in  logical  order  [Cunha  99].  These  techniques  can  be  applied  to  a  wide 
range  of  problems  including  scheduling  and  networking.  “Good”  solutions  to  highly 
complex  problems  can  be  found  within  a  reasonable  amount  of  time.  Methods 
like  genetic  algorithms  lend  themselves  to  parallel  application  due  to  the  need  to 
evaluate  each  member  of  the  population  at  each  generation. 

Advantages  of  heuristic  techniques  lie  in  their  ability  to  locate  solutions 
to  combinatorial  optimization  problems  with  greater  efficiency  than  implicit 
enumerations  techniques.  Moreover,  these  techniques  more  easily  accommodate 
the  discontinuities  and  non-linearities  of  real-world  problems  than  do  gradient- 
based  techniques  [McKinney  94]  [Ritzel  94] .  A  disadvantage  of  these  techniques 
are  the  need  for  large  populations  to  be  evaluated  over  many  generations.  In 
computational  fluid  dynamics,  the  cost  of  analyzing  the  design  for  a  numerical 
simulation  can  be  expensive,  therefore  genetic  algorithms  are  usually  used  with 
less  computationally  expensive,  lower  fidelity  models  [Burg  99]. 

Gradient-based  Method 

Gradient-based  methods  are  conceptually  simple  and  deterministic.  Many 
optimization  algorithms  use  gradients.  Gradient-based  methods  use  both  first 
derivative  ( gradient )  information  and  second  derivative  ( Hessian )  information 
and  are  based  on  the  simple  fact  that  f(x)  increases  or  decreases  in  the  direction  d 
according  as  the  directional  derivative  [VF(x)]'d  is  positive  or  negative.  That  is, 
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the  fastest  way  to  find  an  extremum  is  to  move  along  the  gradient.  Obtaining  these 
derivatives,  however,  can  be  expensive,  and  the  information  may  be  inaccurate. 
There  are  several  gradient-based  methods.  Some  of  the  more  popular  ones  are 
the  method  of  steepest  descent,  conjugate  gradient  methods,  Newton’s  method, 
quasi-Newton  method,  trust  region  models,  and  Gauss-Newton  method. 

The  method  of  steepest  descent  is  often  used  to  calculate  the  next  guess 
xn+l  when  the  only  available  information  is  F(xn )  and  VF(f').  The  primary 
concept  is  that  the  optimal  search  direction  pn  at  xn  is  the  gradient  or  steepest 
descent  direction  VF(ln),  therefore  the  value  of  F(xn )  is  decreased  by  moving 
in  this  direction  with  small  step  sizes.  Drawbacks  to  this  method  are  the  need 
for  one  or  more  additional  function  evaluations  per  design  iteration  and  the  linear 
convergence  rates.  Also,  due  to  poor  convergence,  this  method  may  not  be  efficient 
for  realistic  design  optimization  algorithms  as  many  design  iterations  may  be 
required.  The  gradient  of  the  function  is  defined  as  follows. 

VA(i)«^,l<Kn  (A.l) 

If  ak  denotes  the  optimal  step  length  resulting  from  searching  along  the  direction 
dk  =  — V/(xfe)  starting  from  the  point  xk,  then  the  updated  estimate  of  the  solution 
is  computed  as  follows. 

xk+1  =  xk  -  akVf(xk),  k  >  0  (A. 2) 
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If  the  expression  for  the  gradient  is  available,  then  the  derivative  of  F(a)  needed 
for  the  bisection  method  can  be  computed  using  the  chain  rule  as  follows. 

F'(a)  =  (dk)TVf(xk  +  adk)  (A. 3) 

The  conjugate  gradient  method  uses  F{xn )  and  WF(xn)  like  the  steepest 
descent  method  but  also  has  an  additional  constraint.  The  current  search  direction, 
pn ,  must  be  conjugate,  or  perpendicular,  to  the  previous  design  variable  search 
directions.  For  quadratic  functions,  the  conjugate  gradient  method  identifies 
the  optimal  solution  exactly  within  the  number  of  design  variable  iterations. 
Drawbacks  to  this  method  are  that  optimal  step  size  is  required  to  achieve  efficient 
convergence;  therefore  multiple  function  evaluations  are  necessary  to  determine 
optimal  step  size.  Also,  for  design  optimization  using  high-fidelity  simulation, 
the  conjugate  gradient  method  is  infeasible  due  to  the  high  computational  cost  of 
the  function  evaluation.  The  convergence  rate  scales  with  the  number  of  design 
variables  [Burg  99]. 

Newton’s  method  is  another  gradient-based  method.  It  uses  the  Hessian  to 
update  the  design  variable.  This  method  generally  yields  quadratic  convergence 
once  the  iterates  have  moved  within  the  convergence  region.  A  drawback  to  this 
method  is  the  need  for  the  Hessian  matrix.  This  method  is  rarely  used. 

For  highly  nonlinear  functions,  the  trust  region  model  uses  the  knowledge 
that  the  best  search  direction  as  the  step  size  tends  towards  zero  is  the 
gradient  direction.  The  search  direction  is  approximated  as  a  combination  of  the 
gradient  direction  and  the  Newton  search  direction.  This  method  has  received 
much  attention  as  a  means  to  improve  the  convergence  of  implicit  flow  solvers. 
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Drawbacks  to  this  method  are  the  computational  cost  of  the  function  evaluations 
and  the  need  to  calculate  the  Hessian  matrix. 

The  quasi-Newton  method  gradually  builds  an  approximation  to  the  Hessian 
matrix  based  on  information  at  the  previous  iterations  and  uses  the  Hessian  to 
determine  the  search  direction  and  step  size.  A  requirement  for  this  method  is 
that  the  Hessian  matrix  remain  symmetric  and  positive  definite  and  that  the  new 
Hessian  satisfy  the  secant  equation  Bn+lsn  =  yn.  The  “best”  update  formula 
for  this  method  is  the  Broyden-Fletcher-Goldfarb-Shanno  (BFGS)  according  to 
Broyden  [Broyden] .  These  update  methods  are  often  used  to  achieve  super-linear 
convergence  rates. 

The  Gauss-Newton  method  approximates  the  Hessian  by  assuming  the  second 
order  term  in  V2A  is  negligible.  This  assumption  restricts  the  objective  function 
to  a  least-squares  function.  That  is,  given  a  function 

N 

F(x)  =  r2  (cc)  =  Rt(x)R(x )  (A. 4) 

1=1 

Ti{x)— N  residual  functions 
x—a  vector  of  variables 
R(x)= a  vector  of  residual  functions 

This  method  uses  the  structure  of  the  function  F(x)  to  approximate  the  Hessian 
matrix.  Hence,  it  is  able  to  produce  super-linear  convergence.  The  first  derivative 
of  the  function  is 

8F  4-  d n  dRT  „ 

R—  =  Y1 2H~ri  =  2^j — R 

OXk  OXk  OXk 

1=1 


(A.5) 
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and  the  gradient  is 


YFlx)  =  2  VRtR 


(A.6) 


The  second  derivative  is 


d2F  _  29rf  drj_  +  y^  ^  <92n  _  2dRrdR_  +  2  <92Ar 

dxkXj  dxk  dxj  dxkdxj  1  dxk  dxj  dxkdxj 


(A. 7) 


and  the  Hessian  Matrix  is 


V2F(£)  =  xVRtWR  +  2  S'  (A. 8) 

where  S  is  a  matrix  whose  entries  are  -Jr-^—R.  The  term  S  is  assumed  to  be 

OX^OXj 

relatively  small  compared  to  (Vi?n)2  Vi?"  near  the  local  extrema,  hence  it  can  be 
ignored.  Using  this  assumptions  allows  for  p n  to  be  obtained  via 

2(VRn)TVRnApn  =  -2  (VRn)TRn  (A.9) 

Burg  notes  that  though  the  assumptions  made  do  not  apply  to  his  work,  the 
assumption  that  S  can  be  ignored  appears  to  be  reasonable,  since  the  Gauss- 
Newton  method  obtained  good  design  improvement  results  [Burg  99] . 

Sensitivity  Analysis 

Estimating  the  design  space  gradient  can  be  computationally  expensive  because 
it  would  require  an  additional  steady  state  calculation  for  each  design  variable. 
Using  sensitivity  analysis  can  reduce  this  cost.  The  design  space  gradient  is 
estimated  by  using  information  gained  from  differentiating  the  system  of  governing 
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equations.  There  are  two  basic  types  of  sensitivity  analysis  -  the  discrete  approach 
and  the  continuum  approach.  The  discrete  approach  takes  analytical  derivatives 
of  the  discretized  equations  with  respect  to  the  independent  variables  (Q,  %,  (3). 
The  continuum  approach  calculates  derivatives  directly,  based  upon  the  continuous 
governing  equations,  by  using  the  method  of  material  derivatives.  The  governing 
equations  are  discretized  and  then  differentiated  with  the  discrete  approach.  In 
the  continuous  approach  the  governing  equations  are  differentiated  and  then 
discretized. 

Discrete  sensitivity  analysis  derives  sensitivity  equations  by  differentiating  the 
discretized  system  of  governing  partial  differential  equations  (which  involves  the 
Jacobian).  Many  implicit  simulation  models  use  this  same  Jacobian  matrix,  and 
the  matrix  equation  solution  subroutines  can  be  used  for  the  analysis.  Explicit 
codes  do  not  use  the  Jacobian.  Therefore,  it  must  be  generated  for  sensitivity 
analysis.  The  complex  Taylor’s  series  expansion  (CTSE)  can  be  used  to  generate 
Jacobian  matrix  from  the  discretized  system  of  equations  when  the  Jacobian  does 
not  already  exist.  The  CTSE  can  be  used  to  generate  numerically  exact  Jacobians 
for  the  implicit  codes  as  well. 

The  continuous  approach,  also  known  as  the  continuous  adjoint  approach,  uses 
variational  calculus  to  derive  the  adjoint  equations.  The  objective  function  F, 
the  governing  partial  differential  equations  P,  and  the  boundary  conditions  are 
functions  of  the  flow  variables  Q,  the  computational  domain  or  grid  y,  and  the 
design  variables.  Forming  the  Lagrangian,  taking  the  variation,  and  satisfying  the 
governing  equations  for  a  steady  state  problem  yields 
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SL 


d  L  sr~^  dL 

gf*  +  E  m5l3i 


Sft 


(A.10) 


or 


5L_ 

SPi 


dL  5\  dL 

dx  8  Pi  +  dpi 


dL_ 

dPi 


(A.ll) 


I Q  fixed 

L(Q0),  x(P  +  eAPi),  P  +  eiAPi,  Xj,  rfc) 


-  L(Q(p),x(p),p,xj,rk) 


2A  Pi 


Since  the  partial  differential  equations  and  the  boundary  conditions  are  satisfied 


dF  dL 


Q  fixed 


(A. 12) 


With  the  continuous  approach  a  system  of  partial  differential  equations  is 
differentiated  to  yield  a  system  of  adjoint  equations,  which  is  solved  in  addition 
to  the  governing  equations.  The  difficulties  with  this  method  are  the  need  to 
rewrite  analysis  code  to  solve  a  different  set  of  partial  differential  equations,  the 
computational  cost  of  solving  the  adjoint  equations,  and  the  non-trivial  derivations 
of  the  adjoint  equations.  In  addition,  the  adjoint  variable  equations  must  be  re¬ 
derived  when  changes  are  made  to  the  math  model  [Burg  99] . 

Following  Burg  [Burg  99] ,  the  discrete  approach,  or  direct  differentiation,  solves 
a  discretized  system  of  equations  resulting  from  the  governing  partial  differential 
equations  being  differentiated  at  each  node  in  the  computational  domain  and  a 
linear  matrix  equation.  The  design  space  derivative  of  the  objective  function, 


given  by  F(Q0),  X0),  P),  is 


124 


dF  _  dFT  dQ  dFT  dX  OF 
d(3i  dQ  dPi  +  dX  dPi  +  dpi 


(A.13) 


In  this  equation  dE.  is  the  total  variation  of  F  with  respect  to  the  design  variable 


Pi.  and  can  be  calculated  since  the  dependency  of  F  on  Q,  X,  and  P  is 


OF  9£  onH  df- 
dQ  1  ax’  d/3i 

explicit.  can  be  estimated  by  finite  differencing  the  results  of  the  grid  generation 
code  or  by  differentiating  the  explicit  dependencies  between  the  grid  and  the  design 
variables  with  hand-differentiation,  complex  Taylor  series  expansion,  or  ADIFOR 
[Bischof  92],  Estimation  of  the  vector  ^  with  finite  differences  would  require  an 
additional  steady-state  simulation.  This  term  also  appears  in  the  equation  for  the 
derivative  of  the  discretized  system  of  governing  equations  W (Q (P) ,  X(P) ,  P)  =  0 


or 


dM^_dVVdQ_  dVVfrx  dW_ 
dPi  dQ  dPi  +  dX  dPi  +  dPi 


(A. 14) 


As  in  equation  A.13,  all  terms  except  for  ^  can  be  calculated  without  the  need 
for  a  steady-state  simulation.  Equation  A.  14  along  with  A.13  are  used  to  estimate 
the  design  space  derivative  4^. 

For  the  adjoint  variable  formulation  of  discrete  sensitivity  analysis,  equation 
(A.  14)  is  multiplied  by  an  adjoint  vector  A  and  added  to  equation  (A.13).  This 
yields 


dF_  _  (d_F^_  xTdW\  dQ_  ( _ 

dPi  V  dQ  +  dQ  )  dPi  +  V  dX 


,  dFT  v 

+  -x - 1-  A 


TdW\  dX  dF  TdW 


dX  J  dPi  +  dPi+  ^  dPt 


(A.15) 
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The  coefficient  of  is  set  equal  to  zero  to  avoid  having  to  calculate  for  .  That 

is 


dFT 

~dQ 

dWT 

~W 


+ 

A 


AJ 


■  dW 
~dQ 
d  F 
~dQ 


=  0 


(A. 16) 


A  must  be  calculated  only  once  for  each  objective  function  and  is  independent 
of  the  design  variables  for  the  adjoint  variable  formulation.  The  adjoint  variable 
formulation  scales  with  the  number  of  functions  and  constraints  and  is  appropriate 
for  problems  with  only  one  objective  function  and  a  few  constraints.  The  direct 
formulation  scales  with  the  number  of  design  variables. 

For  the  direct  formulation  of  discrete  sensitivity  analysis,  or  quasi-analytic 
method,  the  vector  ^  is  obtained  by  solving  equation  (A.  14)  directly  or 

dW  dQ  _  _dW  dX  _  dW  _  _dW 

dQ  df3i  dx  d(3i  d(3i  d(3i 

The  vector  ^  must  be  calculated  for  each  design  variable.  This  makes  the 
computational  cost  of  the  direct  formulation  greater  than  that  of  the  adjoint 
variable  formulation  when  the  number  of  design  variables  is  greater  than  the 
number  of  objective  function  and  constraints.  However,  when  using  a  least-squares 
objective  function  and  the  Gauss-Newton  optimization  algorithm  the  number  of 
objective  functions  to  be  differentiated  is  much  larger  than  the  number  of  design 
variables  since  this  algorithm  views  each  term  in  the  objective  function  as  a 
separate  residual  function.  For  this  reason,  the  quasi-analytic  method  is  more 


O  fixed 


(A.17) 


cost  effective. 
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The  advantages  of  the  discrete  sensitivity  analysis  are  that  few  modifications 
are  needed  for  implicit  codes.  The  flow  variables  are  driven  to  steady-state,  and  the 
flow  solver,  the  Jacobian  matrix  dih  and  the  matrix  inverse  routines  are  available 
for  implicit  codes.  Also,  the  computational  cost  for  discrete  sensitivity  analysis  is 
smaller  than  for  the  continuous  approach  due  to  the  need  to  perform  the  pertinent 
calculations  only  at  steady-state.  A  disadvantage  of  the  discrete  sensitivity  analysis 
is  the  need  to  generate  highly  accurate  Jacobian  matrices.  Also,  most  implicit 
flow  solvers  only  approximate  the  Jacobian  matrix  by  neglecting  friction  and 
turbulence  models,  making  assumptions  about  upwinding  methods  and  using  first 
order  schemes  [Piasecki  97].  For  explicit  solution  methods,  the  Jacobian  matrix  is 
not  computed  as  part  of  the  simulation  and  is  unavailable  for  discrete  sensitivity 
analysis.  The  Jacobian  matrix  must  be  computed. 

In  his  work,  Burg  discusses  accuracy  issues  related  to  the  discrete  approach. 
He  notes  that  error  can  be  introduced  into  the  calculations  of  the  design  space 
derivatives  via  several  terms.  The  partial  derivatives  of  the  objective  function 
F,  the  Jacobian  matrix  and  the  vector  must  be  determined.  Also, 

Q  fixed 

the  discrete  approach  assumes  that  the  residual  vector  W (Q)  has  been  successfully 
driven  to  zero  and  that  the  matrix  inversion  operation  is  highly  accurate.  This  may 
not  be  true  if  the  iterative  solvers  are  prematurely  terminated.  Burg’s  work  uses 
hand-differentiation  to  calculate  the  partial  derivative  of  the  objective  function  F 
and  the  complex  Taylor’s  series  expansion  to  calculate  the  Jacobian  matrix.  Burg 
also  notes  that  the  choice  of  grid  and  solution  methods  for  the  continuous  adjoint 
approach  affects  the  accuracy  of  the  adjoint  variable  A.  As  a  result,  the  accuracy 
of  the  design  space  derivative  is  affected.  [Burg  99] 
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